Y F Deng

Prof Yu-Jun J Cui
email: yujun.cui@enpc.fr

A M Tang

X L Li

X Sillen

An Experimental Study on the Secondary Deformation of Boom Clay

Keywords: Boom clay, oedometer test, secondary deformation behavior, mechanically dominated rebounding, physico-chemically dominated swelling

Boom clay formation, a deposit of slightly over-consolidated marine clay that belongs to the Oligocene series in the north east of Belgium, has been selected as a possible host material of nuclear waste disposal. In this context, the long-term deformation behaviour of Boom clay is of crucial importance in the performance assessment of the whole storage system. In this study, low and high pressure oedometer tests are carried out; the e-log σ' v (void ratiologarithm of vertical effective stress) and e-log t (void ratio -logarithm of time) curves obtained are used to determine the compression index C c *, swelling index C s * and secondary deformation coefficient C α during both loading and unloading. The relationship between C α and the effective stress ratio (σ' v /σ' c , vertical effective stress to pre-consolidation stress) is analysed, and it is observed that C α increases linearly with log σ' v / σ' c . Examination of the ratio of C α /C c * for various soils shows that the secondary deformation behaviour of Boom clay is similar to that of shake and mudstone. The relation between C α and C c * is linear; but the relation between C α and C s * is bi-linear. The bi-linearity observed is related to two different mechanisms: the mechanically dominated rebounding and the physico-chemically dominated swelling.

Introduction

Boom clay formation, a thick deposit of slightly over-consolidated marine clay has been selected as a possible host material of nuclear waste disposal in Belgium. In this context, its volume change behaviour, especially its secondary deformation behaviour is essential for the safety of the whole storage system, and therefore needs to be investigated in depth.

The consolidation of fine-grained soils has been commonly described by the primary consolidation and the secondary consolidation. The former refers to the soil volume change due to water pressure dissipation whereas the latter refers to the soil volume change due to the evolution of soil fabric and soil-water interaction. In the past decades, many studies were conducted to correlate the secondary deformation coefficient (C α ) during loading with other soil characteristics. [START_REF] Walker | Secondary settlement in sensitive clays[END_REF] showed that C α varied with the ratio of vertical effective stress (σ' v ) to the pre-consolidation pressure (σ' c ), with the largest C α at a stress slightly higher than σ' c . This was confirmed by other studies on various soils [START_REF] Brook | Secondary compression of soft clay from Ballina[END_REF]Yilmaz and Saglamer, 2004;[START_REF] You | Compression and consolidation behavior under different sampling, loading methods[END_REF][START_REF] Zhu | Consolidation behaviors of soft soil in Nansha[END_REF]Shirako et al., 2006;[START_REF] Suneel | Compressibility characteristics of Korean marine clay[END_REF]Costa and Ioannis, 2009). [START_REF] Walker | The prediction of consolidation rates in a cemented clay[END_REF] found that the secondary deformation coefficient (C α ) during loading has a linear relationship with the compression index (C c ) over the full range of stress applied. This C α -C c relation was further investigated by many other researchers [START_REF] Mesri | The C α /C c concept and Ko during secondary compression[END_REF][START_REF] Mesri | Secondary compression of peat with or without surcharging[END_REF][START_REF] Abdullah | On the physicochemical aspects of compacted clay compressibility[END_REF][START_REF] Al-Shamrani | Application of the C α /C c concept to secondary compression of Sabkha soils[END_REF][START_REF] Brook | Secondary compression of soft clay from Ballina[END_REF][START_REF] You | Compression and consolidation behavior under different sampling, loading methods[END_REF][START_REF] Feng | Consolidation behavior of a soft mud treated with small cement content[END_REF][START_REF] Tan | The behavior of compression and consolidation for clays[END_REF][START_REF] Mesri | Primary compression and secondary compression[END_REF][START_REF] Zhang | Study on secondary consolidation deformation of soft clay[END_REF][START_REF] Zhu | Consolidation behaviors of soft soil in Nansha[END_REF]Costa and Ioannis, 2009;[START_REF] Feng | Experimental study on secondary behavior of soft clays[END_REF][START_REF] Mesri | Compression of granular materials[END_REF] on various soils (intact clays, remoulded clays, clays treated with lime or cement, and sands); the results confirmed the observation by [START_REF] Walker | The prediction of consolidation rates in a cemented clay[END_REF]. [START_REF] Mesri | Settlement of embankment on soft clays[END_REF] defined four groups of soils according to the value of the ratio C α /C c (Table 1). Some other correlations were also attempted between C α /C c (or C α ) and soil physical properties such as the liquid limit w L , plastic limit w P and plasticity index I p [START_REF] You | Compression and consolidation behavior under different sampling, loading methods[END_REF][START_REF] Suneel | Compressibility characteristics of Korean marine clay[END_REF].

Although the secondary consolidation behaviour of soils has been widely investigated, there have been few studies on the stiff Boom clay, especially on the unloading path that represents the situation of the soil in vicinity of excavated galleries. In the present work, consolidation tests are performed in both low and high pressure oedometers on Boom clay taken from the sites of Essen and Mol, Belgium. Loading and unloading are run in steps and the secondary deformation coefficient C α is determined for each step. Furthermore, the relations between C α , the effective stress ratio (σ' v /σ' c ), compression and swelling indexes (C c * and C s *) are analyzed. The main objective is to study the variations of C α during unloading and the mechanisms involved in these variations. Note that the use of high pressure oedometer in this study allows studying the variation of C α at large stress ratio σ' v /σ' c , indispensable for deeply located soil as Boom clay (223 m deep in Mol and about 240 m in Essen). Moreover, the introduction of parameter C c * and C s * allows analysing the soil compression behaviour with a non-linear loading-unloading curve. Note also that, to the authors' knowledge, there have been no studies before focusing on the variations of C α during unloading.

Soil studied

The soil studied was taken by coring at the sites of Essen and Mol, Belgium. The location of the two sites is shown in Figure 1 [START_REF] De Craen | Geochemical analyses of Boom Clay pore water and underlying aquifers in the Essen-1 borehole[END_REF]. The Essen site is situated in the north east of Belgium, about 60 km far from the underground research laboratory (URL) at the Mol site. After being taken from the borehole, the cores were sealed in plastic tubes having ends closed and transported to the laboratory. Five soil cores of 1-m length and 100-mm in diameter from Essen and one soil core of 0.5-m length and 100-mm in diameter from Mol were studied. The details of these cores are shown in Table 2, with the corresponding depth, member, unit mass of solids (ρ s ), liquid limit (w L ), plastic limit (w P ), plasticity index (I p ), water content (w 0 ) and void ratio (e 0 ). There are three cores taken from the Putte member (Mol, Ess75 and Ess83) and three cores from the Terhagen members (Ess96, Ess104 and Ess112). The geotechnical identification parameters of the cores from Essen are similar: ρ s = 2.64 -2.68; w L = 62 -78%; w P = 25 -33%; I P = 36 -45. The values are also close for both water content and void ratio: w 0 = 27.2 -29.7, e 0 = 0.700 -0.785. For the core from Mol, the values of ρ s, w L , w P and I P are similar to that of the cores from Essen, but the water content and void ratio are lower than the cores from Essen, showing that Boom clay from Mol is denser.

Experimental techniques

Both low pressure (0.05 -3.2 MPa) and high pressure (0.125 -32MPa) oedometer tests were carried out following the French standards [START_REF] Afnor | Sols : reconnaissance et essais: essai de gonflement à l'oedomètre, détermination des déformations par chargement de plusieurs éprouvettes[END_REF][START_REF] Afnor | Geotechnical investigating and testing, Laboratory testing of soils, Part 5: Incremental loading odometer test[END_REF] on the six Boom clay cores.

The tests in low pressure oedometer aim at studying the loading-unloading behavior of the soil near the excavation gallery, whereas the tests in high pressure oedometer aim at studying the compression behavior in large stress level (far from the excavation gallery). The soil samples were prepared by trimming and had 50-mm in diameter and 20-mm in height. In the following, high pressure oedometer test is named Oedo1 while low pressure test is named Oedo2. Note that the high pressure oedometer has the same principle as the standard low pressure oedometer; the main difference is that in high pressure oedometer two amplification levels were used with a ratio of 1:10 for the first level and 1:5 for the second level (see Figure 2). In other words, the frame of high pressure oedometer allows multiplying the applied weight by 50, which leads to a maximum force of 12 tons. In the experiments, the minimum and maximum applied weights are 5 N and 1280 N, leading to a minimum and maximum vertical pressure of 0.125 MPa and 32 MPa respectively for a sample of 50 mm diameter.

The soil specimen was installed in the cell with dry porous stones. Prior to circulation of the synthetic water which has the same chemistry as the in-situ pore water [START_REF] Cui | Investigating the time dependent behavior of Boom clay under thermomechanical loading[END_REF] in the drainage system, a confining pressure equal to the estimated in situ stress was applied. This prevents the soil swelling during re-saturation which may modify the soil microstructure and as a result the soil mechanical properties [START_REF] Delage | Suction effects in deep Boom clay samples[END_REF].

The in situ stress of the soil was estimated using Eq. 1:

0 ' 0 u h v - = γ σ (1)
where σ' v0 is the in situ effective vertical stress; γ is the mean unit weight of the soil above the depth considered, taken equal to 20 kN/m 3 following the data of De Craen et al. (2006); h is the depth of the soil core (see Table 2); u 0 is the in situ pore pressure estimated from the ground water level that is assumed to be at the ground surface. The ' 0 v σ values determined for Ess75,Ess83,Ess96,Ess104,Ess112 and Mol are 2.20,2.27,2.40,2.48,2.56 and 2.23 MPa,respectively. For a reason of convenience, σ' v0 in both low pressure and high pressure oedometers was set at 2.40 MPa for all tests.

Experimental results

Compressibility behavior

Figure 3 presents the loading-unloading-reloading stages and the corresponding changes in vertical displacement in test Ess75Oedo1. Before the re-saturation phase, a loading from 0.125 to 2.4 MPa was applied to reach the in situ stress state. The soil sample was then re-saturated using synthetic water. The subsequent unloading-reloading stages were as follows:

unloading from point A (2.4 MPa) to B (0.125 MPa); loading to C (16 MPa); unloading to D (0.125 MPa); loading to E (32 MPa) and unloading to F (0.125 MPa). Common results were obtained in terms of vertical displacements, i.e. compression upon loading and rebounding upon unloading. Note that the French standards -AFNOR (1995,2005) were applied as regards the deformation stabilization for all odometer tests: stabilization is achieved when the displacement rate is lower than 0.01mm/h. The low pressure oedometer test Ess75Oedo2 shows a quasi elastic behavior with narrow unloading-loading loops. A deeper examination shows, however, that the reloading curve from II to III is not linear in the plane e-log σ' v . This non-linearity can be also observed on the curve of test Ess75Oedo1 on the reloading paths from B to C and from D to E. Note that the results from the tests on other cores are similar to that shown in Figure 4. Obviously, it is difficult to determine the pre-yield stress σ' y using the Casagrande method on these curves. In addition, this pre-yield stress, if any, does not correspond to the pre-consolidation pressure σ' c (σ' y is much lower than σ' c ): σ' c is equal to 2.4 MPa at point A, but when reloading from B to C, σ' y seems to be much lower, about 1 MPa. For this reason, in the following analysis only For the determination of secondary deformation coefficient C α , standard method is used based on the e -∆log t plot. The determination of C c *, C s * and C α is illustrated in Figure 5.

Note that C α =∆e/∆log t is negative when loading and positive when unloading.

Relation between C

α α α α and σ σ σ σ' v /σ σ σ σ' c
Figure 6 shows the variation of C α versus the stress ratio σ' v /σ' c for all the six cores, identified by both low pressure and high pressure oedometer tests. It appears that during loading stage C α ranges mostly from 0 to 0.01 especially when the vertical effective stress is lower than the pre-consolidation stress. Some points beyond 0.01 can be observed when the stress ratio is greater than 2. For core Ess112, C α is very small and close to zero when the stress ratio is less than 1. For the unloading stages, C α ranges mostly from 0 to -0.01 for all tests when the stress ratio is greater than 0.1. On the contrary, when the stress ratio is less than 0.1, C α is less than -0.01. These observations lead to conclude that more significant secondary consolidation takes place at higher stress ratios (σ' v /σ' c > 1) upon loading and more significant secondary swelling takes place at lower stress ratios (σ' v /σ' c < 0.1).

Except the results of core Ess112 during loading, all other results show that C α increases almost linearly with the stress ratio in the semi-logarithmic plane, for both loading and unloading stages. This is different from the results reported by other authors [START_REF] Walker | Secondary settlement in sensitive clays[END_REF][START_REF] Brook | Secondary compression of soft clay from Ballina[END_REF][START_REF] Yilmaz | Evaluation of secondary compressibility of a soft clay[END_REF][START_REF] You | Compression and consolidation behavior under different sampling, loading methods[END_REF]Zhu et is greatly acknowledged for its financial support. The first author is grateful to the National Science Foundation of China for its support (No 50908049). Table 1. Soil classification according to the values C α /C c [START_REF] Mesri | Settlement of embankment on soft clays[END_REF] Table 2. Geotechnical properties of the soil cores studied

List of Tables

List of Figures

Figure 4

 4 Figure 4 presents the compression curve (void ratio versus logσ' v ) of test Ess75Oedo1, together with the compression curve of test Ess75Oedo2 in low pressure odometer. In test Ess75Oedo2, after the re-saturation using synthetic water under 2.4 MPa stress, unloading was performed from point I (2.4 MPa) to point II (0.05 MPa), then loading to point III (3.2 MPa) and finally unloading to point IV (0.05 MPa).

  σ' c is used and its determination is based on the stress history: σ' c is the maximum stress applied in the oedometer tests. For instance, σ' c = 16 MPa for the paths C->D and D->E; σ' c = 32 MPa for the path E->F; σ' c = 3.2 MPa for the path III->IV is 3.2 MPa. Since the e-log σ' v curves of Boom clay are not linear during unloading and unloading, especially for the low pressure oedometer tests, it is difficult to use an unique compression index (C c ) and swell index (C s ) to describe the compression and swelling behavior. Hence the above two indexes are determined stage by stage, and renamed C c * and C s *, respectively. It should be pointed out that if the e-log σ' v curves for the three stages (i.e. before pre-yielding, after pre-yielding and unloading) are linear, the C c * becomes the same as C c and C s * becomes the same as C s .

Figure 7

 7 Figure 7 also shows that during unloading, a bi-linear relation between C α and C s * can be observed: the turning point at a C s * value around 0.1. This turning point can be considered as an indicator of changes from mechanical dominance to physical-chemical dominance in terms of volume changes: when C s * is less than the value at the turning point, the clay shows a mechanically dominated rebounding; by contrast, when C s * is larger than the value at the turning point, the clay shows a physico-chemically dominated swelling. This particular behaviour during unloading was also observed in other works: Delage et al. (2007) and Le et al. (2011) conducted compression tests on unsaturated Boom clay with suction monitoring, and observed that during unloading the soil suction increased slowly in the beginning and then rapidly when the vertical stress decreased down to a threshold value; Cui et al.(2002)

Figure 1 .Figure 5 .

 15 Figure 1. Locations of the sampling sites (De Craen et al., 2006) Figure 2 Sketch of high pressure oedometer Figure 3. Vertical effective stress and displacement versus elapsed time (Ess75Oedo1) Figure 4. Compression curves from oedometer tests (Ess75Oedo1 and Ess75Oedo2) Figure 5. Determination of parameters C c *, C s * and C α Figure 6. C α versus stress ratio σ' v /σ' c . (a) Core Ess75; (b) Core Ess83; (c) Core Ess96; (d) Core Ess104; (e) Core 112; (f) Core Mol Figure 7. C α versus C c * and C s *. (a) Core Ess75; (b) Core Ess83; (c) Core Ess96; (d) Core Ess104; (e) Core 112; (f) Core Mol Figure 8. C α versus C c * Figure 9. C α versus C s *

Figure 1 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 9 .

 1345679 Figure 1. Locations of the sampling sites (De Craen et al., 2006)

  Based on the classification criterion in Table1, one can conclude that Boom clay falls in the zone of shake and mudstone whose C α /C c * value ranges from 0.02 to 0.04.

[START_REF] Afnor | Geotechnical investigating and testing, Laboratory testing of soils, Part 5: Incremental loading odometer test[END_REF][START_REF] Shriako | Secondary compression behavior in standard consolidation tests[END_REF][START_REF] Suneel | Compressibility characteristics of Korean marine clay[END_REF] Costa and Ioannis, 2009)

who observed that the relation between C α and log σ' v /σ' c during loading is rather convex.

4.3. Relation between C c * or C s * and C α α α α

Figure

7

shows the variations of C α with C c * and C s * for all the cores. It appears that C α increases linearly with C c *, with a slope ranging from 0.019 to 0.029. Moreover, this linear relation is independent of the state of consolidation: for a given core, all the points below and beyond σ' c are on the same line.

Mesri et al.

(1994)

analyzed the secondary consolidation behaviour of many soils, and gave the correlation between the secondary consolidation coefficient (C α ) and the compression index C c as shown in Table

1

. From the results obtained on Boom clay, it appears that the ratio C α /C c * falls in a narrow range from 0.019 to 0.029. In order to have a mean value, all the results during loading are gathered in Figure

8

, in terms of variations of C α versus C c *. A value of 0.024 is identified for the ratio C α /C c *.

Table 1 .

 1 Soil classification according to the values C α /C c[START_REF] Mesri | Settlement of embankment on soft clays[END_REF]

	Material	C α /C c
	Granular soils including rockfill	C α /C c = 0.02 ± 0.01
	Shake and mudstone	C α /C c = 0.03 ± 0.01
	Inorganic clays and silts	C α /C c = 0.04 ± 0.01
	Organic clays and silts	C α /C c = 0.05 ± 0.01
	Peat and muskeg	C α /C c = 0.06 ± 0.01

Table 2 .

 2 Geotechnical properties of the soil cores studied

	344
	Core	Depth (m)	Member	ρ s (Mg/m 3 )	w L (%)	w p (%)	I p (%)	w 0 (%)	e 0
	Ess75	218.91-219.91 Putte	2.65	78 33 45 29.7 0.785
	Ess83	226.65-227.65 Putte	2.64	70 33 37 27.2 0.730
	Ess96	239.62-240.62 Terhagen	2.68	69 33 36 26.5 0.715
	Ess104 247.90-248.91 Terhagen	2.68	68 29 39 27.7 0.700
	Ess112 255.92-256.93 Terhagen	2.67	62 25 37 27.3 0.755
	Mol	223	Putte	2.67	68 26 42 23.6 0.625

Acknowledgements

ONDRAF/NIRAS (The Belgian Agency for Radioactive Waste and Enriched Fissile Materials)


Pierre Bénard

Jingwan Lu

Forrester Cole

Adam Finkelstein

Joëlle Thollot

Active Strokes: Coherent Line Stylization for Animated 3D Models

Keywords: NPR, line drawings, temporal coherence, snakes

Figure 1: Stylized animation of a galloping horse. From left to right: line samples are extracted from a 3D model; active strokes track the samples; brush paths are attched to the strokes and stylized as circular arcs; two more frames of animation exhibit temporal coherence.

Introduction

Line drawings created by artists often include stylistic effects such as textured brushes and graceful, curved strokes. These effects reflect the motion of an artist's tool, and as such are fundamentally 2D in appearance. The main challenge in recreating such effects for computer animation is to create lines that smoothly track the 3D scene while maintaining the desired 2D texture and shape. This paper presents a method to create such a set of temporally coherent, stylized lines from an animated 3D model. We term these lines active strokes, because they adapt and extend active contours ("snakes") [START_REF] Kass | Snakes: Active contour models[END_REF]] to stylized rendering. Active strokes automatically update their topology and shape as the underlying 3D scene animates, and provide coherent parameterization for texture and shape abstraction (Figure 1). Our approach performs at interactive framerates, making it suitable for artistic exploration, games, interactive visualizations and illustrations.

Stylization effects such as texture and shape abstraction require that lines be represented as connected curves. Suitable curves can sometimes be extracted and stylized by traversing the 3D scene in object space [START_REF] Northrup | Artistic silhouettes: a hybrid approach[END_REF][START_REF] Grabli | Programmable rendering of line drawing from 3D scenes[END_REF]. In order to produce smooth, coherent animation, several approaches have been proposed to preserve the parameterization of object-space lines over time [START_REF] Kalnins | Coherent stylized silhouettes[END_REF]Bénard et al. 2010;[START_REF] Buchholz | Spatio-temporal analysis for parameterizing animated lines[END_REF][START_REF] Karsch | Snaxels on a plane[END_REF]. However, for models of moderate complexity like the Stanford bunny shown in Figure 2, objectspace line extraction may produce many noisy, short segments that confound stylization with spatial and temporal coherence. Moreover, object space lines lack natural level of detail (LoD) control, and offer no straightforward mechanism to depart from the underlying surface in support of 2D shape abstraction effects.

Active strokes avoid these issues by operating in image space, using line samples output by filter-based line extractors (e.g. [START_REF] Saito | Comprehensible rendering of 3-D shapes[END_REF][START_REF] Raskar | Image precision silhouette edges[END_REF][START_REF] Lee | Line drawings via abstracted shading[END_REF]). Since no natural connectivity exists between these samples, active strokes first vectorize the line samples using snakes. Like traditional snakes, our snakes conform to the changing shape of the features in the scene. Unlike traditional snakes, they also need to adapt their topology and overall arrangement during animation. As the snakes split and merge to track the line samples, their parameterization may change discontinuously. Active strokes therefore carry independently parameterized brush paths on top of the snakes to maintain consistent stretches of parameterization from one frame to the next (Figure 3). Since the snake handles tracking of motion in the scene, the brush paths may deviate substantially from the position of the underlying line to support shape abstraction.

The key innovation of our method is the introduction of an independent 2D entity, the active stroke, that approximates lines extracted from the 3D model and persists from frame to frame. We present new mechanisms to add, remove, trim, extend, split and merge these strokes to match animated feature lines. Our method provides level of detail control, operates on line samples from any type of line extractor, and supports new stylization effects such as shape abstraction for lines. We demonstrate our method with a line rendering system that achieves interactive frame rates and provides a broad range of styles.

The definitive version is available at http://diglib.eg.org/

Related Work

Line stylization. Image processing is the simplest, most efficient way to extract lines from 3D models. The seminal work of [START_REF] Saito | Comprehensible rendering of 3-D shapes[END_REF] filters depth and normal maps to extract contours and creases. It produces line drawings with natural coherence and LOD in 2D, but does not support much stylization or abstraction due to the raster nature of the output. Generalizing the method of [START_REF] Lee | Line drawings via abstracted shading[END_REF], [START_REF] Vergne | Implicit brushes for line-based rendering[END_REF] combine 2D local fitting with spatially-varying convolution to produce more complex styles. But their implicit definition is without line parameterization, preventing the application of stroke texture or abstraction. Our approach supports image-based line extraction, but recovers connectivity and parameterization in support of such effects.

Object space methods extract connected paths on the surface of the model. These paths have obvious parameterization for each individual frame [START_REF] Grabli | Programmable rendering of line drawing from 3D scenes[END_REF]], but the parameterization can change abruptly during animation due to changes in line topology. [START_REF] Kalnins | Coherent stylized silhouettes[END_REF] proposed the first complete method to enforce coherent parameterization by using the connectivity of smooth silhouettes [START_REF] Hertzmann | Illustrating smooth surfaces[END_REF] to propagate stylization. This approach works well for simple models with few lines, where the mapping between object and image space is almost one-to-one. Models of moderate complexity like the Stanford bunny, however, produce many, tiny line fragments that collide in image space, thwarting effective propagation (Figure 2). The method of Bénard et al. [2010] enforces a common parameterization for nearby lines (thereby making short segments behave as parts of a longer line) but this solution is generally suitable only for nearby lines that are also parallel. The method of [START_REF] Karsch | Snaxels on a plane[END_REF] tracks snake-like snaxels on the 3D object rather than re-extracting at each frame, avoiding propagation of parameterization. Because the snaxels live in object-space, however, they cannot easily provide LoD control or abstraction, and do not appear to address the aforementioned line fragmentation problem.

Another contribution of Kalnins et al. is the use of multiple brush paths per silhouette. The parameterization of each brush path is independently optimized to find a compromise between uniformity in screen-space and in object space. Our brush paths behave similarly, with a few differences. Their brush paths both track 3D motion and transmit parameterization from frame to frame. Our snakes handle the former role, so that our brush paths can be abstracted, deviating far from the snakes without compromising the quality of tracking (Figure 3). In addition, we employ the self-similar line artmap (SLAM) textures of Bénard et al. [2010] to allow extra flexibility in the brush path parameterization.

To date most line stylization approaches operate in an online fashion, i.e., without knowledge of the future positions of the lines. [START_REF] Buchholz | Spatio-temporal analysis for parameterizing animated lines[END_REF] present a method for parameterizing lines across an entire animation, allowing lines to anticipate topology changes and merge smoothly. While presented in the context of object-space occluding contours, thus sharing the difficulties of other object-space approaches, the spatio-temporal concept is independent and could be applied to our active strokes, though we have not explored that option here.

Active contours in NPR. Introduced by [START_REF] Kass | Snakes: Active contour models[END_REF], snakes are widely used in computer vision and medical imaging applications because of their robustness to noise. They have also been applied in several NPR contexts: line extraction [START_REF] Karsch | Snaxels on a plane[END_REF], painterly rendering [START_REF] Hertzmann | Paint by relaxation[END_REF]], and video stylization [START_REF] Agarwala | Snaketoonz : A semi-automatic approach to creating cel animation from video[END_REF][START_REF] Agarwala | Keyframe-based tracking for rotoscoping and animation[END_REF]]. As mentioned, the snaxels approach of Karsch and Hart works in object space, where the topology of the contours generators forms a set of closed loops. On the other hand, methods that use snakes to delineate image space region boundaries need to deal with a network of open curves, and thereby take on the need to deal with more complex topology changes over time. The systems of Agarwala et al. segment and track color regions in videos with B-spline snakes for coherent stylization. The user manually draws the initial outlines of the contours and the outlines are propagated to track features in the video. These systems were designed for processing video offline and can therefore rely on heavy-weight optimization and manual correction. In contrast, our approach is designed to track motion field in animated 3D scenes where the camera path is not known in advance, for example games. Therefore it must be interactive, fully automatic, and cannot rely on "future" frames. On the other hand, we benefit from the quality of our input motion -the motion field extracted from 3D is generally less noisy than optical flow from video.

Figure 3: Each brush path (dotted) is defined relative to a snake (solid), but can deviate substantially to support abstraction effects as shown in Figure 1. As the knot turns, the snake topology evolves; for example the orange snake merges with the purple snake (top). Coherent arc-length parameterization is preserved by maintaining the beige brush path over the purple snake. Line drawing LOD. Line level of detail (LoD) control is necessary to produce clean, abstract drawings from complex models. Image-space approaches provide a natural limit on line complexity (image pixels), but it is often desirable to further reduce density. Object-space lines allow arbitrary line density, but provide opportunities for simplification both during stylization [START_REF] Grabli | Density measure for line-drawing simplification[END_REF][START_REF] Cole | Directing gaze in 3d models with stylized focus[END_REF]] and as a post-process [START_REF] Barla | Geometric clustering for line drawing simplification[END_REF][START_REF] Shesh | Efficient and dynamic simplification of line drawings[END_REF]. A common approach is to extract lines, measure their density in screen space, then either omit or merge lines in overly dense areas. [START_REF] Grabli | Density measure for line-drawing simplification[END_REF] and [START_REF] Barla | Geometric clustering for line drawing simplification[END_REF] found that good results could be achieved by removing only overlapping lines that also align in orientation. We adopt a similar, simple definition of the line coverage that allows efficient processing and explicit control over spatial and directional proximity. Finally, by generating brush paths that are independent of the feature lines, we can easily merge lines rather than simply omitting them, unlike approaches such as [START_REF] Cole | Directing gaze in 3d models with stylized focus[END_REF].

Overview

The key idea of our approach is create an independent set of curves (active strokes) that track and approximate unorganized feature line samples and provide a coherent parameterization for stylization. Active strokes use snakes in order to track and vectorize the input feature samples. On top of the snakes an active stroke carries a set of brush paths that propagate coherent parameterization and give the final shape of the strokes. Figure 4 illustrates our pipeline.

Our method takes as input a set of feature samples with 2D position and 2D tangent, and for animated scenes, 2D velocity. These samples can be extracted by any method that provides an orientation at each pixel (e.g. [START_REF] Lee | Line drawings via abstracted shading[END_REF][START_REF] Vergne | Implicit brushes for line-based rendering[END_REF]), but the smoother and more coherent the input samples, the better the result. We have found that lines extracted with steerable quadrature pair filters [START_REF] Freeman | The design and use of steerable filters[END_REF] provide especially good input. We describe our application of these filters in Appendix A.

For the first frame only, we construct a set of snakes conforming to these feature samples via a stochastic process described in Section 5.

For successive frames, the snakes track and vectorize the input samples. To obtain a good tracking, we maintain two properties for our snakes: coherence, meaning snakes should evolve continuously from frame to frame as they track the features of the model; and accuracy, meaning snakes should faithfully represent the shape of the underlying linear features extracted from the scene. For vectorization, we maintain the following properties: coverage, each feature line should be covered by one snake and each snake should cover a feature; simplicity, snakes should maintain simple, smooth and clean topology; length, snakes should be as long as possible while respecting the previous criteria.

To achieve coherence, snakes from the previous frame are advected to follow the motion field in the scene (Section 4.1). Next the snakes are relaxed to conform to the features in the current frame (Section 4.2) in support of accuracy. After advection and relaxation, the goals of coverage, simplicity and length will generally have been violated. The ideal solution would involve a complex non-linear optimization, difficult to achieve in a system designed for interactive frame rates. Instead, we adopt a greedy approach and update each snake individually, using heuristics to move towards a better solution (Section 5). We obtain efficient vectorization by beginning with a good guess adapted from the previous frame.

Finally, we build on top of the snakes a set of brush paths that provide a coherent parameterization across frames (Section 6.2). By relying on the tracking behavior of the snakes, the brush paths can deviate widely from the underlying feature lines to achieve stylization and abstraction goals (Section 6.1). The brush paths are then textured and drawn to produce the final illustration.

Tracking Feature Lines

The first step of our approach is to track the feature lines available in each frame of the animation using snakes: 2D polylines that act like elastic strands, simultaneously attempting to satisfy external attraction forces and internal smoothing forces. In our approach the snakes are attracted to feature lines and track them from frame to frame. At each frame, we move the snakes according to the motion field in the scene (Section 4.1) and then relax them to conform to the feature lines in the current frame (Section 4.2).

Advection

At initialization, each snake vertex is associated with a 2D feature sample (Section 5). During advection, the snake vertex is moved according to the projected 3D motion field of the scene. The 3D motion field can be computed either by back projecting the 2D samples in the camera at frame i -1, reprojecting these 3D positions into the camera at frame i and computing the difference or by using an optical flow computation when the geometry is not available.

After advection, lines from frame i -1 are nearby (but not perfectly aligned with) features in frame i. At this stage, the method of [START_REF] Kalnins | Coherent stylized silhouettes[END_REF] performs an explicit search to associate line samples with features in frame i. Instead, we use relaxation to attach the advected snakes to features in frame i, as described next.

Relaxation

As described by [START_REF] Kass | Snakes: Active contour models[END_REF], snakes are parametric 2D curves v(s) = (x(s), y(s)) , s ∈ [0, 1] that minimize the energy E, defined as the sum of smoothing and data terms Eint and Eext:

E = 1 0 Eint(v(s)) + Eext(v(s)) ds
The internal energy Eint ensures continuity (first-order membrane term weighted by α) and smoothness (expressed as second-order thin-plate term weighted by β):

Eint(v(s)) = 1 2 (α(s)|v ′ (s)| 2 + β(s)|v ′′ (s)| 2 )
with v ′ and v ′′ the first and second derivatives. In practice, we use constant α and β parameters (both equal to 0.001 in our experiments).

Addressing the accuracy goal, the external energy Eext attracts the snake to areas of interest -in our case, feature lines that are projected and rasterized into a binary image I. This image is filtered by a Gaussian kernel whose gradient is the attraction field:

Eext(v(s)) = -| ∇ [Gσ(v(s)) ⋆ I(v(s))] | 2
The size σ of the filter kernel controls the width of the attraction field around the contour. A larger filter makes relaxation more robust to errors in advection, at the price of accuracy in tracking. In our experiments we use σ = 8px for screen-resolution images.

Following the approach of Kass et al. we optimize the energy E by gradient descent. The continuous snake v(s) is discretized as n vertices vi and the energy is iteratively minimized by a semi-implicit Euler scheme. At each frame, a matrix of internal forces is constructed and applied iteratively, updating the contours positions and their associated external forces at each iteration. This scheme is simple and fast so long as the advection step provides a good starting condition. But when advection fails to produce a good initial guess it is likely due to extreme motion of the model. In this case, temporal coherence may not be necessary or even expected.

To maintain good sampling, especially while zooming, the snake resolution is updated after several iterations of relaxation. Via the approach of [START_REF] Delingette | New algorithms for controlling active contours shape and topology[END_REF] we ensure an almost uniform sampling of the curves, with a maximum distance of δmax = 6px.

Snakes have a natural tendency to shrink like a rubber-band due to their internal forces. This shrinking tends to produce inaccurate results where snakes round off corners or fail to reach line endpoints. To overcome this behavior, we apply spring forces along each snake edge so that they maintain almost constant length if no tangential external force is applied. The spring forces are implemented as an extra term in the external energy, similar to the approach of [START_REF] Mcinerney | Topologically adaptable snakes[END_REF]:

fi→i+1 = k (Li -vi+1 -vi ) vi+1 -vi vi+1 -vi
and the symmetric case for fi→i-1 where k is the stiffness parameter of the spring (k = 1 by default) and Li is the rest-length equal to vi+1 -vi at the previous frame.

Vectorizing Feature Lines

During animation snakes must be created, deleted, or modified to preserve a good approximation of the underlying feature lines. Constructing snakes from the feature samples is equivalent to vectorizing line art. This is a difficult problem even for still images, and our scenario also requires that the vectorization be temporally coherent. Fortunately, we can exploit the smooth tracking behavior of the snakes to ease the vectorization problem. Snakes from the previous frame that have been advected and relaxed provide an excellent starting point for vectorizing the samples for the current frame, because the features they track are coherent in image space.

However, starting with these snakes, we still need to address two challenges. The first challenge is that after advection and relaxation we generally have not met our unit coverage goal. Several snakes may now overlap in image space where a single snake would suffice. Also, new lines may have appeared in regions where no snake existed before, or conversely a snake may remain where feature lines no longer exist. The second challenge is snake topology. Unlike our scenario, typical formulations for parametric snakes need not address the issue of automatically inferring their connectivity, usually by relying on user input. Starting from this initial connectivity, some approaches have been proposed to modify topology of the contours during segmentation tasks [START_REF] Mcinerney | T-snakes: Topology adaptive snakes[END_REF][START_REF] Delingette | New algorithms for controlling active contours shape and topology[END_REF][START_REF] Ji | Robust topology-adaptive snakes for image segmentation[END_REF]. We take inspiration from these methods to define a set of heuristic rules that allow the snakes to change their topology and overall arrangement during animation so as to match the underlying feature lines.

We organize these heuristics as six operators -delete, split, trim, extend, merge, insert -which are applied sequentially in a greedy fashion on the vertices of the snakes, until they can no longer modify the snakes. The order in which the operators are applied is important. Deletion and splitting occur first, so as to allow subsequent trimming to produce clean junctions. Extension also benefits from following splitting, and must precede merging to allow snakes to join across a gap. Insertion of new snakes occurs to address features without coverage after all the other operations have had the opportunity to adjust the snakes (especially extension). So in summary, we delete as much as possible, trim as much as possible, and so on; and finally insert new snakes where extra coverage is needed.

On the other hand, the order in which the snakes are processed for each operation has relatively less influence on the quality of the final result. Nevertheless, we process snakes in order from longest to shortest, except for the trim operator for which we use the opposite order, with the goal of maximizing the length of the longer snakes.

Operations

Our goal is to cover each feature line with exactly one snake. Recall that each feature line is actually composed of an unorganized set of feature samples S that are typically sparse in the image. We define coverage mean that a feature sample si with tangent #» Ti has a snake vertex vj nearby whose local orientation #» tj is similar:

covered(si) ≡ ∃vj : si -vj < rcover | #» Ti • #» tj | > θc
The constants rcover and θc control, respectively, the apparent line density and the tendency for snakes to match the orientation of the underlying features. A related notion is footprint(vj) which describes the set of feature points covered by vj:

footprint(vj) = si ∈ S | ( si -vj < rcover) ∧ (| #» Ti • #» tj | > θc)
Initialization. For the first frame of an animation we start by creating new snakes as follows. We randomly pick an uncovered feature sample from which we grow a snake as far as possible using the extension operator described below. This process is repeated until it can no longer improve coverage. Though not providing optimality guarantees, this approach is simple and fast, and it offers a reasonable first guess for subsequent processing.

Delete. Snake vertices that are not covering any feature sample are deleted. This mechanism may remove an entire snake where a feature has become occluded, or may cut away portions of a snake in case of partial occlusion. 5a) To accurately track sharp features of a 3D model, a snake is split at vertex vi if a sharp angle appears there, found by:

Split. (Figure

c split ≡ # » ti-1 • # » ti+1 < θ ⊥
where # » ti-1 and # » ti+1 are the tangents at two neighboring vertices and θ ⊥ is a threshold constant for how much the snake is allowed to bend at a vertex. Trim. (Figure 6a) Having deleted portions of snakes that do not cover features, and split at sharp corners, we next trim snakes in regions where parallel lines provide more coverage than needed. Several strategies might be possible, for example attempting to merge lines where they nearly overlap. However, we have found a simple policy that trims lines at their endpoints to be effective for achieving unit coverage. Specifically, the three criteria for removing the endpoint vi from a snake are that it is too close to a vertex vj of a different snake, they are close to parallel, and that by removing vi we will not leave a feature uncovered that would otherwise be covered by vi:

c trim ≡    vi -vj < r trim | #» ti • #» tj | > θ footprint(vi) ⊆ ∪ k =i footprint(v k )
where #» ti and #» tj are their respective tangents, and r trim and θ are threshold constants for how close and parallel these parts of the snakes may become before trimming occurs. 6b) Next we extend snakes at every endpoint where doing so will offer coverage for otherwise uncovered feature points whose tangents align locally with the snake. Specifically, beyond the current endpoint vi we will consider appending a new endpoint vj at the location of every uncovered feature sample, based on four criteria:

Extend. (Figure

c extend ≡        ¬ covered(vj) vj -vi < r extend #» ti • #» tj > θ | #» tj • # » Tj| > θ
where #» ti and #» tj are the tangents at the snake vertices, # » Tj is the tangent of the feature sample, and r extend a search radius.

In general there may be many candidate feature locations for which the extension criteria above are satisfied. Therefore, we compute a goodness score gj that describes how well each candidate extends the snake and then choose the candidate with the highest score. The score favors aligning the tangents and a relatively short extension, as follows: where the two terms in parentheses have been normalized to the range [0, 1] and α is a parameter that balances between them. 5b) To simplify the snake topology, two snakes are merged if they are locally parallel and their endpoints are close:

g extend = α min( #» ti • #» tj , | #» tj • # » Tj |) -θ 1 -θ +(1-α) r extend -vj -vi r extend -rcover

Merge. (Figure

cmerge ≡    #» τi • #» tj > θ #» τj • #» ti > θ vj -vi < rmerge
These criteria may lead to multiple conflicting candidate mergers, for instance where three endpoints meet in a Y shape. Therefore, in a fashion similar to the extension process described above, we assign a goodness score that favors straightness and closeness:

gmerge = α min( #» τi • #» tj , #» τj • #» ti ) -θ 1 -θ + (1 -α) rmerge -vj -vi rmerge -rcover
where vi and vj are the endpoints of two distinct contours, #» ti and #» tj their tangents, and #» τi and #» τj are tangents of hypothetical merger segments, shown in Figure 5b. We perform as many non-conflicting mergers as possible, in order of descending score.

Insert. This final operation works like the initialization step described above, except that it occurs after all the other operations. Its role is to insert new snakes where new features have appeared that have not been addressed by, for example, extending existing snakes.

Practical considerations

Threshold Constants. The scope of action of the four operations that trim, extend, split and merge the snakes are given by threshold constants for various radii and angles. Several relationships between them are important. In order for the snakes to establish proper coverage between samples: rcover > δmax/2. In order to permit extension to occur: r extend > rcover. In order to enforce hysteresis between split and merge events: θ ⊥ < θ . Other constants trade off between accuracy and performance. In our experiments we have found the values effective for the full range of results shown: rcover = 10px; r trim = 5px; r extend = 20px; rmerge = 20px; θc = θ = 0.85; θ ⊥ = 0.5; α = 0.8. Acceleration structure. During advection and coverage updates, we often need to access the neighbors of snakes vertices or feature samples in the coverage radius. To expedite local search, we use the data structure of [START_REF] Delingette | New algorithms for controlling active contours shape and topology[END_REF]: a randomized sparse regular grid (implemented by a hash table) containing all the snakes vertices and feature samples. With a grid cell size of 2 rcover, we enumerate all points in a disc of radius rcover by walking four cells.

Brush paths shape and parameterization

Snakes approximate and track the underlying feature lines from frame to frame. Based on this temporally coherent layer, we build a set of brush paths that are in charge of (1) providing the shape and position of the final strokes, and (2) of computing a temporally coherent parameterization.

Brush path geometry

We represent a brush path as a 2D polyline with vertices associated 1 -1 to a sequence of vertices on a single snake (Figure 3). The position of brush path vertices are described relative to their associated snake vertices to take advantage of the snake coherence during motion. The user controls a target length lmax (potentially infinite) for the brush path allowing for short brush paths if desired. In addition, brush paths can have "overshoot" that extends their end segments.

Smoothing in time

During rapid camera movements, feature lines may appear or disappear abruptly. As snakes persist across frames, we can smooth out these events by assigning a lifespan to the snake vertices, as follows. Each snake vertex has a lifespan κ, initially null. At each frame, if a vertex persists, its lifespan increases: κ = max(κ + 1, κmax) with κmax a constant defined by the user (typically 3 our experiments). Conversely, if a vertex disappears, lifespan decreases (κ = κ -1). This lifespan is transmitted to the brush path vertices and can be mapped to different style attributes, such as line opacity or thickness. That way, brush strokes fade in and out progressively during κmax frames.

Shape Abstraction

The shape of the brush path can deviate substantially from the shape of its corresponding snake, which allows highly abstract styles. As a proof of concept we show two extreme examples that never been animated before: straight brush paths and circular arcs (Figure 7).

Drawing with straight brush paths implies to partition each snake in a set of approximately linear segments. We compute this partition using dynamic programming, similarly to the approach of [START_REF] Mccrae | Sketching piecewise clothoid curves[END_REF] for fitting sketched clothoid curves.

Between two frames, the initially straight brush paths are propagated following the motion of the underlying snakes. However the shape of the snakes is likely to have changed. Consequently, the brush paths may not be straight anymore, and we may even need to break a straight path into two or more segments to better match the new shape. Such a decision is taken using the same dynamic programming approach applied independently on each brush path. In addition, wherever the policy described in Section 6.2 might consider merging brush paths, we ensure that the dynamic programming solution would not immediately split the merged path.

Conveniently, the abstraction algorithm for straight brush paths is easily adapted for circular arcs. We use the "modified least squares" method of [START_REF] Umbach | A few methods for fitting circles to data[END_REF] to determine a best arc and compute the fitting error for a sequence of vertices.

Slowing motion When using strong abstraction, temporal artifacts may still occur. For example a circular arc can suddenly switch from concave to convex, or a single segment can break into two. To soften these transitions, we provide the following optional smoothing mechanism. Suppose the smoothed position of a path vertex at frame i is si. At frame i + 1, rather than sending the vertex to its target position ti+1, we move it a step in that direction d = (ti+1 -si):

si+1 = si + min(1, dmax/ #» d ) #» d
where dmax specifies a maximum step size. Note that si and ti are represented as offsets relative to the corresponding snake vertex and not absolute positions. Thus our mechanism leaves unchanged the motion of the model while smoothing the abstracted paths. The smoothing is applied on the endpoints of segments, whereas every arc vertices are processed (meaning that the latter shape deviates from exact circles during smoothed transitions).

Brush path parameterization

We need to parameterize these paths with temporal coherence during animation. Snakes track features smoothly evolving in 3D. In absence of visibility events, their intrinsic parameterization correlates with the 3D motion of the feature lines. But it does not correspond to screen-space arc-length, which induces compression artifacts. Moreover, occlusions and disocclusions produce sliding and popping. To avoid such effects, we recompute the parameterization at each frame based on the previous one, and we use the self-similar line artmap (SLAM) textures of Bénard et al. [2010] to avoid compression when zooming. To handle topological eventsmerging, typically -we allow multiple brush paths per snake, and forbid merging of two brush paths with different parameterization.

Parameterization fitting.

Similarly to the 2D coherence strategy of [START_REF] Kalnins | Coherent stylized silhouettes[END_REF], we assume that a brush path represents the interaction between a drawing tool and 2D paper. Thus it should be uniform in screen-space, implying that its parameterization T is linear with respect to arc-length s: T (s) = ρs + φ. Besides, T should evolve according to the motion of the underlying snake. To account for these two goals, we follow the approach of Kalnins et al., a linear least squares fit of T , based on the parameter values from the previous frame propagated to the brush path of this frame (via the associated snake). However, our approach offers an extra degree of freedom: SLAM textures allow us to handle gracefully variations in phase and slope while preventing sliding artifacts, especially when zooming. Also note that our method does not need a fitting algorithm robust to outliers, such as RANSAC used by Bénard et al.

[2010], because we are not mixing the parameterizations of multiple brush paths.

Topological events. Brush paths rely on the connectivity of their associated snakes. When a snake is updated so are its brush paths. For example, snake resampling (Section 4.2) resampling of the paths. When a snake is split (Section 5.1), the brush paths covering the break are also split; likewise for trimming, extension and deletion. However, when two snakes merge, their brush paths are not merged automatically, as their endpoints may differ in position and parameterization. Instead, a new merging operator ensures that the paths are continuous spatially and in parameter value. If not, both paths are kept (Figure 3). Otherwise, the two brush paths are joined and their mutual parameterization becomes:

T (s) = ρs + φ1 where ρ = ρ2l2 + φ2 -φ1 l1 + l2 To avoid brush path fragmentation over time, we also propose a leveling mechanism, similar in spirit to the "healing" process of [START_REF] Kalnins | Coherent stylized silhouettes[END_REF]. It progressively pushes the parameterization of two reasonably close brush paths to their common mean T , encouraging a fusion in subsequent frames. Note that we can explicitly deal with the periodicity of the parameterization during this process, which increases its efficiency.

Results

Snakes and their accompanying brush paths provide a robust, temporally coherent and evenly parameterized set of curves in screen space. Thanks to these qualities, we can explore styles that previous approaches were unable to apply on animated 3D scenes of moderate complexity (Figure 7 and Figure 11).

Style examples

We can produce standard drawing styles by using SLAM textures together with stroke attributes such as color and thickness. Tapering [START_REF] Northrup | Artistic silhouettes: a hybrid approach[END_REF] is provided and can be computed per brush path or per snake.

Offsets We can also apply displacement mapping to the vertices of each brush path, exploiting the fact that brush path geometry is persistent across frames. This mechanism supports, for example, the loopy effect added to the knot or the woman in Figure 7. [START_REF] Kalnins | WYSIWYG NPR: drawing strokes directly on 3D models[END_REF] screen space offsets using Markov random fields and apply them on the brush paths geometry. In contrast, we synthesize a texture of offsets using the SLAM synthesis approach of Bénard et al. [2010] which provides similar details at varying scales.

Overdraw The unit coverage computation at the snakes level not only prevents line clutter but also supports a form of "overdrawing," i.e. the overlap of a controlled number of brush paths on the same snake. The user specifies a target overdraw Ωt ∈ R > 1. We measure the current overdraw Ωc between two snake vertices as the number of brush paths overlapping this edge. For every brush path endpoint, we stochastically decide whether to extend, trim, or leave it unchanged, based on the probabilities:

P (extend) = max(0, γ extend (1 -Ωc/Ωt)) P (trim) = max(0, γ trim (1 -Ωt/Ωc))
where γ extend and γ trim weight how rapidly the overdraw moves towards its target during animation. We can create even sketchier styles by combining overdraw with a small target length and random offsets per brush paths.

Speed lines Brush paths are persistent across frames and have their own linear parameterization. We take advantage of these two properties to produce speed lines (Figure 8) inspired by lines appearing in hand drawn comics.

Speed lines are created as new brush paths formed by chaining together copies of seed vertices over a series of frames. These seeds are chosen at locations undergoing noticeable motion, i.e. regions of brush paths where the velocity #» v has a strong component opposite to the projected surface normal. At every brush path vertex we compute V = -#» v • n. We also compute the mean µV and standard deviation σV of V all the vertices. A vertex is defined as a seed when V > t σV + µV with t a threshold that controls how much of the model receives these effects (we use 1 2 ).

Speed lines can then be rendered as any other brush path using the effects described earlier. We generally find that it helps to trim the end of the speed lines after a given length.

Performance

Our implementation is primarily CPU-based, except for the 2D line extraction, blur for relaxation, and final rendering which rely on the GPU. The system performs at interactive rates on a recent machine (Intel Core 2 at 2.40 Ghz with 4 GB of RAM and a GeForce GTX 260) for models of moderate complexity (ranging from 20 to 30 FPS on the models shown in this paper). A major bottleneck is the readback from the GPU to CPU, which is performed at several stages. Techniques such as stream compaction [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF]] could be used to reduce this cost. The sparse matrix inversion required for the relaxation step is performed efficiently using the CHOLMOD library [START_REF] Chen | Algorithm 887: Cholmod, supernodal sparse cholesky factorization and update/downdate[END_REF].

Comparison with Previous Approaches

Our method is the first to allow a coherent parameterization of image-space feature lines. Implicit Brushes [START_REF] Vergne | Implicit brushes for line-based rendering[END_REF] was the best available image space approach for line stylization.

We obtain more complex styles with parameterization based effects (regular textures, tappering) and shape abstraction.

Previous object-space methods [START_REF] Kalnins | Coherent stylized silhouettes[END_REF]Bénard et al. 2010;[START_REF] Buchholz | Spatio-temporal analysis for parameterizing animated lines[END_REF][START_REF] Karsch | Snaxels on a plane[END_REF] rely on visible stretches of object space feature lines as input, and their brush paths can not be longer than these features. The impact of brush path length on the final rendering can be seen in Figure 9. The short, unpredictable paths of the silhouette segments are clearly visible, while the snakes provide long brush paths that correspond with the major shape features of the model.

An alternative solution to our snake based approach could be to chain the small pieces of line together with an involved chaining algorithm (e. g., [START_REF] Grabli | Programmable rendering of line drawing from 3D scenes[END_REF]) before computing coherence using [START_REF] Kalnins | Coherent stylized silhouettes[END_REF] or [START_REF] Buchholz | Spatio-temporal analysis for parameterizing animated lines[END_REF]. However, achieving interactive performances is very unlikely with such an approach. Moreover, vectorizing fully disconnected image space samples wouldn't be possible.

As a nice side effect, working in image space allows us to provide a LoD mechanism. While zooming in and out, the constant coverage radius ensures a unit line density for a given range of orientations (Figure 10). That way, brush paths appear and disappear naturally according to the camera motion as can be seen in the accompanying video.

Figure 10: Zoom out on an octopus depicted with object space lines (left) and active strokes (right). The line screen space density varies with our approach avoiding line clutters (bottom).

Limitations

Our method deals with image-space lines. We therefore have little knowledge of the 3D informations of the scene. It mostly has an impact on occlusions and disocclusions. For these visibility events a snake does not know what are the best feature samples to follow. It creates a tearing effect where during 2 or 3 frames the snake zigzags between two feature lines. This is the price to pay to have an imagebased approach allowing for efficient frame-rates that would not have been possible to obtain with an object-space approach.

Another limitation of our approach is that for efficiency, decisions about snakes and brush paths are made locally. There is no explicit mid-or high-level understanding of the scene. For example, tracking T-junctions might be a way to improve the behavior at occlusions. As a result, our approach is well-suited for feature lines that are view-dependent and extracted from organic shapes, rather than fixed features such as those found in architectural scenes.

Conclusions and future work

This paper describes an approach for coherent evolution of feature line rendering for animated 3D scenes. Feature lines are approximated by snakes that advect with the motion of the scene. The snakes adjust their arrangement and topology using a new set of operators, in order to match the changing scene. Attached to the snakes, brush paths carry parameterization from one frame to the next with coherence, and can therefore be stylized and abstracted with temporal coherence. We introduce several new stylization effects shown in Section 7 and in the accompanying video.

Our pipeline is generic enough to be applied on video, by using vision algorithms to extract an approximated version of the data that we need: optical flow, feature lines, and depth. This endeavor would be particularly timely with the advent of commodity cameras that supply depth and optical flow along with color information.

Finally, while this system cannot optimize globally over an entire animation sequence (since we assume the frames are not known in advance) it might be possible to anticipate one or two frames in the future using something like a Kalman filter on the camera path and then use this information to broaden the scope of the optimization. A Line Extraction with Steerable Filters [START_REF] Freeman | The design and use of steerable filters[END_REF] describe an image contour extractor using steerable quadrature pair filters. The quadrature pair responds equally to step-edges and fine lines in the image, providing a consistent, intuitively appealing set of line samples. The steerable formulation provides smooth, precise estimation of orientation, which we use for the tangent at each contour pixel.

We apply the contour extraction to a shaded image, allowing us to capture lighting as well as shape effects (similar to [START_REF] Lee | Line drawings via abstracted shading[END_REF]). The image is optionally augmented with depth information to ensure that depth discontinuities are captured by lines. The input to the method is a scalar image I(x, y) = L(x, y) + βz where x and y are the screen space pixel coordinates, L(x, y) is any shaded image, and β is a parameter that controls the influence of depth z. We apply to I the second-order quadrature filter pair contour detector of [START_REF] Freeman | The design and use of steerable filters[END_REF]. Finally, we perform nonmaxima suppression using the estimated orientation to find the set of line samples. The parameters of the method are the filter kernel size, a noise threshold, and β.

The filter pair can be computed very efficiently on the GPU using the steerable formulation with seven basis filters. Each basis filter can be expressed as the composition of a Gaussian and a small derivative kernel, where the Gaussian can be computed with a separable formulation. We use this method for all the results with the exception of Figure 11.

Figure 2 :

 2 Figure 2: Silhouette coherence. Frame 1: (a) Silhouettes in object space appear to be continuous in the image, but (b) are extracted as many disconnected parts. (c) Snakes recover 2D connectivity. Frame 2: (a) The silhouette is visually coherent with frame 1, but is (b) composed of a different set of disconnected parts. Nevertheless, the snakes (c) provide coherence.

Figure 4 :

 4 Figure 4: Stages of our line drawing pipeline, an overview of which is given in Section 3 below.

Figure 5 :

 5 Figure 5: Split and merge operations modify snake connectivity. (a) When a snake bends too sharply around a corner it is split. (b) When two snakes abut smoothly they are joined.

Figure 6 :

 6 Figure 6: Modifying snakes to improve coverage. (a) The endpoint of a snake is trimmed if it is too close and parallel to another snake. When multiple snakes overlap in screen space this operation moves towards unit coverage. (b) The end of a snake is extended when doing so adds coverage to otherwise uncovered features. In concert, the trim and extend operations tend to clean up T-junctions.

Figure 7 :

 7 Figure 7: Stylization examples, from left to right: monkey with circular arcs (above) and fuzzy SLAM texture over toon shading (below); knot with construction lines (above) and overdrawn wiggles (below); woman in two poses and three styles: arcs, loopy offsets, and overdrawn.

Figure 8 :

 8 Figure 8: Speed lines stylized with arrows and dots textures.

Figure 9 :

 9 Figure9: The Jot system of[START_REF] Kalnins | Coherent stylized silhouettes[END_REF] (a) generates many short strokes because the visible silhouettes produce short brush paths. Our approach (b) recovers connectivity using snakes, leading to fewer, longer strokes.

Figure 11 :

 11 Figure 11: Animated models stylized with active strokes: running elephant rendered with circular acrs (left) and segments (middle); animated hand in three poses with SLAM texture and offsets.

Acknowledgements

Thanks to Szymon Rusinkiewicz, Fredo Durand, Sylvain Paris, Pascal Barla and Aaron Hertzman for their insightful comments. Many of the images shown in this paper are rendered from 3D models courtesy of the Princeton Shape Benchmark, the Stanford Computer Graphics Laboratory, the Utah 3D Animation Repository, and the Computer Graphics Group at MIT. This work was sponsored in part by the ExploraDoc program of the region Rhône-Alpes, France.


Christyves Chevallier
email: christyves.chevallier@supelec.fr

Nicolas Fressengeas

Frédéric Genty

Joel Jacquet

Chevallier C • Fressengeas

N • Genty

F • Jacquet

Chevallier C • Genty

Fressengeas N Lmops

Robust design of Si/Si 3 N 4 high contrast grating mirror for mid-infrared VCSEL application

Keywords: High contrast grating mirror, mid-infrared VCSEL, Robust design

teaching and research institutions in France or abroad, or from public or private research centers.

Robust design of Si/Si3N4 high contrast grating mirror

for mid-infrared VCSEL application

Introduction

High contrast grating (HCG) mirrors with a grating period smaller or in the range of the reflected wavelength can be designed to exhibit very high performances with energy only in the 0-order of diffraction. HCGs which are made of a high index grating on top of a low index sublayer can achieve a reflectivity as high as 99.9 % and bandwidths larger than 100 nm together with a strong polarization selectivity between TM and TE mode [1,2]. These characteristics are very interesting for replacing the top Bragg mirror in VCSEL devices [3]. Moreover, the circular VCSEL structure usually presents no polarisation selectivity and the use of a grating should improve the stability of the emit-ted beam [4]. Also, the wavelength of emission of mid-IR electrically-pumped VCSEL is currently limited close to 2.6 µm since the necessary semiconductor Bragg mirrors become as thick as 11 µm with more than 20 quarter wavelength pairs [5,6] making the epitaxial growth control critical. With a reduction of a factor of 5 of the thickness of the top mirror and fabrication tolerances as large as 20 % [7,8], the use of HCGs should allow VCSEL operation at larger wavelengths. Moreover, the HCG dimensions are scalable with the wavelength in the limit of the refractive index dispersion [1] which makes the designs presented in this work scalable for other wavelengths.

In this work, we design a Si/Si 3 N 4 HCG for a VCSEL application at 2.65 µm. First, the optical properties that VCSEL mirrors must exhibit are defined and quantitatively expressed through a quality factor. The dimensions of the HCG structure are then adjusted by an optimization algorithm to maximize the quality factor and meet VCSEL requirements. Finally, the HCG robustness is enhanced through a precise study of the tolerance on the physical dimensions. The final design has thus an optimized reflectivity with several percent of tolerance as an error margin that the technological manufacturing process must respect.

Mirror design

In a previous works [8], Si/SiO 2 structures were designed for applications at λ =2.65 µm. However in the 2.65 -3 µm wavelength range, OH-radicals present an absorption band [9] and could decrease the HCG performances depending on their concentration in the SiO 2 layer [8]. The Si 3 N 4 should be well suitable for the low index material since it exhibits no absorption at these wavelengths [9]. The mirror structure presented in Fig. 1 is made of a silicium grating (n = 3.435) on top of a Si 3 N 4 layer as the low index sublayer (n = 1.99). A Si/Si 3 N 4 quarter wavelength pair has been added below the HCG in order to increase the width of the reflectivity bandwidth [2]. The mirror structure has thus several dimension lengths to define : the etched length L e , the filled length L f , the grating thickness T g and the sublayer thickness T L . Contrary to Bragg mirrors, HCG design rules to achieve high reflectivity are not straightforward and require an optimization process to match the expected properties.

The VCSEL application requires specific properties on the optical performances of the mirror in order to create a cavity of sufficient quality for laser emission. Due to the thin gain region of the VCSEL structure, the reflectivity of the mirror is typically chosen to be as high as 99.5 % for the largest possible bandwidth. However, a 99.9 % reflectivity threshold has been chosen to ensure a fabrication security margin. Moreover, due to their one dimensional symmetry, HCG can be shaped to obtain a polarization selectivity. Thus, the reflectivity threshold of 99.9 % is applied only for transverse magnetic (TM) mode while keeping a transverse electric (TE) reflectivity below 90 %, so that the laser effect should be possible only in TM polarization. The mirror is simulated by the rigorous coupled wave analysis method (RCWA) [10,[START_REF] Rathgen | mrcwa 20080820[END_REF] which computes reflection spectra for TM and TE modes for perfectly squared grating with an infinite number of periods. Calculations were done by using constant index values.

Through a numerical analysis of the computed reflection spectra, the quality of a design is quantitatively evaluated by the use of a quality factor Q [2]:

Q = ∆λ λ 0 1 N λ2 λ=λ1 R T M (λ)g(λ) (1)
This quality factor is written to take into account the previously defined VCSEL requirements with a polarization selectivity performed by selecting the range of wavelengths ∆λ around λ 0 for which R T M is higher than 99.9 % and R T E is lower than 90 %. This selection starts at the aimed wavelength λ 0 and is enlarged both for larger and smaller wavelengths around λ 0 while the reflection conditions are met giving two boundaries : λ 1 and λ 2 . The centering of the bandwidth around λ 0 = 2.65 µm is made by computing a gaussian weighted average of the R T M values on the N points of the bandwidth [λ 1 , λ 2 ]. Finally the influence of the width of the stopband is expressed in this quality factor by multiplying the latter average by the bandwidth itself ∆λ = |λ 2 -λ 1 | and normalized by λ 0 .

Optimization

The optimization of the mirror performance is made by adjusting the grating parameters to maximize the quality factor Q. However, with its numerous local maxima, its non linearity and unknown derivative, this Q factor can only be optimized using a global maximization method. The Differential Evolution algorithm [START_REF] Kroshko | [END_REF]13] has been used as optimization algorithm to find automatically the best set of parameters which maximizes the mirror quality.

The Si/Si 3 N 4 -based structure presented in Fig. 1 has been optimized under constraints by using large and shallow grooves (L f , L e > 500 nm and an aspect ratio T g /L e < 1.1). These limitations are relatively severe for the optimization algorithm since the Si/SiO 2 HCGs are typically defined with narrower patterns of less than 260 nm and large aspect ratio of more than 2.6 [1]. Since HCGs require a vertical etching profile with squared pattern to reach the 99.9 % high reflectivity, the technological constraints were chosen to obtain an experimental square grating profile as close as possible to the theoretical one. However, these values can be easily adjusted before the optimization by the manufacturers according to the etching process used and the machine specifications.

The resulting optimum grating parameters (Table 1) are L e = 829 nm, L f = 522 nm, T g = 899 nm, T L = 896 nm. This mirror design exhibits a 218 nm large bandwidth with R T M > 99.9% together with a high polarization selectivity by keeping R T E < 90% and meets all previously defined VCSEL requirements.

In order to evaluate the feasibility of the design, a tolerance study has been performed on the dimensions of the grating. The tolerance of one parameter is defined as the allowed variation range of this parameter for which the mirror is performant enough for a VCSEL application with R T M > 99.9 % and R T E < 90 %. From a fabrication point of view, the most critical parameters of the grating are the etched length L e and the filled length L f which result from the etching of the slab. The grating thickness can be controlled with more precision during the growth and a selective etching process technique should allow more control on the depth of the etching.

The tolerance measurement of the optimized design has shown large variation range of ∆L e = 68 nm and ∆L f = 24 nm. However the optimum point is located at a boundary of these variation ranges and the allowed error is of only 2 nm on L e and even 0 nm on L f (Fig. 2).

Robustness enhancement

In order to increase the tolerance of the design, a more centered point within the variation range should be used. However, during the fabrication process, errors can be made simultaneously on several parameters at a time, especially on L e and L f if the period was not well defined during the photolithography step for instance. To complete the tolerance study, the variation ranges of the period Λ = L e + L f and the fill factor F F = L f /Λ have been evaluated. These extremum values have been plotted on a 2D map in Figure 3. Thus, the variation range of L e , L f , the period Λ and the fill factor F F = L f /Λ indicates the area of the (L e , L f ) couples where the previously defined VCSEL requirements are met. The optimum point (829, 522) with the largest bandwidth found by Fig. 2 Evolution of the normalized bandwidth versus the etched length Le and the filled length L f . The arrows show the optimum values for these parameters which are both at the upper limit of the variation range.

Fig. 3 The tolerance map represents the variation range for the (Le, L f ) couples while keeping Tg and T L at their optimum value. The points (+) define a polygon (in grey) which vertices correspond to the extrema of the variation range of (Le, L f ) couples where R T M is below 99.9 % or R T E above 90% at λ 0 (∆λ = 0 nm). The center of the largest incircle of the polygon (⋆) indicates the most tolerant point with respect to the fabrication errors on Le and L f . the optimization algorithm is located at the edge of the tolerance area. So, in order to enhance the robustness of the grating with respect to the error of fabrication which could be made on L e and L f , a non optimum point located at the centre of the tolerance area has been chosen with L e = 784 nm and L f = 537 nm.

The resulting HCG, which characteristics are summarized in Table 2, exhibits a 186 nm large bandwidth which is 32 nm less than the optimum design reported in Table 1, but is more robust with tolerance values of ±4 % for L e , ±5 % for L f and ±2 % for T g . These values are not as high as the Si/SiO 2based HCG performances which can exhibit 250 nm large bandwidth with up to 10 % of tolerance. These better performances are due to the higher index contrast of the Si/SiO 2 system (n SiO2 = 1.509).

Conclusion

In this work, a Si/Si 3 N 4 high contrast grating mirror has been designed and exhibit a 186 nm large 99.9 % high reflectivity bandwidth. The dimension of the mirror structure has been selected by an optimization algorithm which maximizes a quality factor specially defined for a VCSEL application. Then, the tolerances of this optimum design have been evaluated and the robustness with respect to the error of fabrication have been improved by adjusting the grating parameters. This mirror, with a thickness of less than 2 µm, grooves larger than 500 nm and several percent of tolerance on the grating parameter, should allow laser operation of VCSEL devices at 2.65 µm.

Fig. 1

 1 Fig. 1 Scheme of the HCG optimized for λ 0 = 2.65 µm exhibiting a 218 nm large bandwidth with Le = 829 nm, L f = 522 nm, Tg = 899 nm, T L = 896 nm.

Table 1

 1 Dimensions and tolerances of the optimum mirror.

	Parameter	Optimum	Minimum Maximum Tolerance
	Le	829 nm	763 nm	831 nm	± 0.2 %
	L f	522 nm	498 nm	522 nm	± 0 %
	Tg	899 nm	897 nm	941 nm	± 0.2 %
	T L	896 nm	806 nm	976 nm	± 8.9 %
	∆λ	218 nm
	∆λ/λ 0	8.3 %

Table 2

 2 Dimensions and tolerances of the robust mirror.

	Parameter	Value	Tolerance	Minimum Maximum
	Le	784 nm ± 4 %	735 nm	814 nm
	L f	537 nm ± 5 %	512 nm	598 nm
	Tg	899 nm ± 2 %	882 nm	923 nm
	T L	896 nm ± 37 %	566 nm	1496 nm
	∆λ	186 nm
	∆λ/λ 0	7 %

Acknowledgment

The authors thank the French ANR for financial support in the framework of Marsupilami project (ANR-09-BLAN-0166-03) and IES and LAAS (France), partners of LMOPS/Supélec in this project. This work was also partly funded by the InterCell grant (http://intercell.metz.supelec.fr) by INRIA and Région Lorraine (CPER2007).