Olivier Gourmel

Adrien Bernhardt

A Gradient-Based Implicit Blend Preprint

Keywords: Descriptors: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling Implicit Curves & Surfaces, Interactive Techniques Implicit surfaces, Blending, Constructive solid geometry (CSG), Geometric modeling, Surface-solid object representations •

We introduce a new family of binary composition operators that solves four major problems of constructive implicit modeling: suppressing bulges when two shapes merge, avoiding unwanted blending at a distance, ensuring that the resulting shape keeps the topology of the union, and enabling sharp details to be added without being blown up. The key idea is that field functions should not only be combined based on their values, but also on their gradients. We implement this idea through a family of C ∞ composition operators evaluated on the GPU for efficiency, and illustrate it by applications to constructive modeling and animation.

INTRODUCTION

Functionally-based implicit surfaces [Bloomenthal 1997b] were introduced in the eighties as promising alternatives to mesh-based and parametric representations. They opened the way to constructive modeling systems where different object components could be smoothly blended, and enabled the production of animations where fluid-like objects could change topology over time through successive merging and splitting. These effects were achieved by combining field functions using a composition operator. Surprisingly, after being enhanced in the nineties by the introduction of convolution surfaces [Bloomenthal 1997b] and of general construction trees [START_REF] Wyvill | Extending the csg treewarping, blending and boolean operations in an implicit surface modeling system[END_REF], little advance was made on such operators in the last decade, leaving four major problems unsolved:

(1) Bulging problem: Implicit blending creates unwanted bulging.

For instance, a shape generated by blending together some implicit cylinders defined in the same plane will necessarily be thicker where the cylinders join, as shown in Figure 1(a-left) where the 'A' is composed of three blended cylinders. The expected result is the 'A' on the right.

(2) Locality problem: Implicit models typically blend at a distance. This is a major problem in modeling applications, where preventing the blend between at least some parts of the models (e.g. between the hand and the head of the character in Figure 1(b)) is necessary. The issue is serious in animation too, where pieces of soft material should be allowed to deform each other and eventually blend, but only where they collide.

(3) Absorption problem: Sharp and thin details are smoothed and they blow up when blended into larger implicit primitives, since they are totally included into the support region of the latter. This problem prevents the creation of thin details such as the eyelashes in Figure 1 (b). This is why implicit surfaces have the reputation of only producing simple, blobby shapes.

O. Gourmel et al. (4) Topology problem: Lastly, although implicit modeling is an excellent way of generating arbitrary topologies, blending often produces material that will fill a hole that a designer intended.

For example the center of the 'A' in Figure 1(a-left). The user has no easy way to control, or at least to predict, the resulting topology.

Solving the above problems is essential for both constructive modeling and animation. While individual solutions have been proposed, none of them solves all these four issues at once. Doing so efficiently opens the way to many applications, from sketch-based, interactive design of 3D shapes (Figure 1(b)) to the interactive animation of deformable solids (Figures 1(c) and 1(d)). We review the state of the art and the contributions before presenting the required technical background and our general solution to these four problems.

Previous work

Implicit surfaces are well known for their ability to combine shapes by simply composing the associated field functions. Composition is generally expressed as a binary operation, except for the n-ary operators; max, which results in a union and sum, which generates smooth blending between the input shapes [START_REF] Sabin | Function representation of solids reconstructed from scattered surface points and contours[END_REF][START_REF] Ricci | A Constructive Geometry for Computer Graphics[END_REF][START_REF] Blinn | A generalization of algebraic surface drawing[END_REF]]. Improving composition operators has been a major topic of research, and modeling techniques were enhanced by the introduction of clean union operators generating an exact union of the input shapes but with a smooth field everywhere else and of blending operators with local shape control [START_REF] Pasko | Function representation in geometric modeling: concepts, implementation and applications[END_REF][START_REF] Barthe | Controllable binary csg operators for "soft objects[END_REF]]. Other operators were developed for animation applications, such as the generation of a contact surface surrounded by bulges when soft objects collide [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF]].

The elimination of the bulging problem was only tackled in two seminal papers: [START_REF] Rockwood | The displacement method for implicit blending surfaces in solid models[END_REF]] proposed an operator, the superelliptic blend, where the range of the blend is parameterized by the cosine of the angle between the field gradients. With this formulation, the blend smoothly transitions to a union where the composed objects' surfaces become collinear. Even though suppressing the unwanted bulges, this operator was C 0 only and increased both the locality problem and the topology problem, as reported in [Bloomenthal 1997a]. The latter discussed alternative solutions for the n-ary composition of skeleton-based primitives: the kernels or the skeletons of convolution surfaces were used to parameterize a switch from sum to max operators. In addition to not being generally applicable, the lack of smoothness of these early solutions limited their usability in the case of successive blends. Our work builds on [START_REF] Rockwood | The displacement method for implicit blending surfaces in solid models[END_REF]] since our new composition operators are smooth functions of the angle between field gradients.

Although identified early [Bloomenthal 1997b], the locality problem was only solved recently in a general setting [START_REF] Pasko | Bounded blending for Function-Based shape modeling[END_REF][START_REF] Bernhardt | Implicit blending revisited[END_REF]. The proposed solutions are based on an extra implicit primitive, the blending volume, whose field is used to interpolate between blending and clean union (such as those of Figure 2). [START_REF] Bernhardt | Implicit blending revisited[END_REF]] automatically generates blending volumes around intersection curves and sets parameters such that small primitives do not blow up, providing a solution to the absorption problem as well, however this method has some drawbacks:

(1) only G 1 continuous field functions are produced leading to curvature artifacts illustrated with reflection lines in Figure 3(c).

(2) Computations involve sampling the intersection curve between the two input surfaces or extracting the curve from the intersection of the meshes used for display. This makes the solution costly and not scalable, especially in applications where many intersection curves would need to be generated at different resolutions, such as in a detailed model.

The new method presented here is a more general, scalable approach and is able to restrict blending to regions where the input shapes intersect. At the same time the topological genus of the resulting shape can be constrained to remain that of the union (topology problem), which none of the previous solutions could achieve.

Among recent uses of implicit surfaces, we can notice several experimental sketch based modeling systems [START_REF] Pasko | Function representation in geometric modeling: concepts, implementation and applications[END_REF][START_REF] Karpenko | Free-form sketching with varational implicit surfaces[END_REF][START_REF] Alexe | Shape modelling by sketching using convolution surfaces[END_REF][START_REF] Wyvill | Sketch based construction and rendering of implicit models[END_REF][START_REF] Bernhardt | Matisse : Painting 2d regions for modeling free-form shapes[END_REF][START_REF] Brazil | Sketching variational hermite-rbf implicits[END_REF] for which our new composition method would be especially useful.

Contributions

We present a generic family of composition operators that brings a unified solution to the four major problems mentioned above. The key idea is to parametrize the blending by the angle between the gradients of the combined field functions, as suggested by [START_REF] Rockwood | The displacement method for implicit blending surfaces in solid models[END_REF]. This enables us to interpolate between a union and a blend in the same composition. Using these operators we generate smooth blending between some parts of input shapes and union elsewhere. Blends near intersections are localized without requiring any extra implicit primitive, which makes the method efficient and scalable. In contrast with early work the operators described using our system are C ∞ continuous, which is a sufficient property to handle successive compositions.

Our blending technique preserves sharp details and gives enough control on the blend localization so that in general, the topology of the composed object can remain that of the union, which makes composition predictable and well suited to constructive modeling frameworks. Our operators can also be tuned to animate progressive separation of soft material while preventing disjoint parts from blending at a distance. They can also generate contact surfaces surrounded by bulges, possibly followed by progressive merging, when soft objects collide.

The method is generic among functionally-based implicit surfaces.

Although it can be applied to discrete fields as well, where f is defined as some interpolation of sample values stored in a grid, we do not address level-set modeling frameworks since the composed field functions we output are not governed by differential equations.

This paper derives our solution for local implicit primitives, such as soft-objects and meta-balls. As explained in the next section, these surfaces are the most challenging to handle since the composed field has to maintain compact support. However, our framework also applies to globally supported primitives, as we illustrate in Section 5.1 with convolution surfaces. Section 2 presents the necessary background and we refer to [Bloomenthal 1997b] for a full description of implicit models.

BACKGROUND

Implicit surfaces: A functionally-based implicit surface S is the C-iso-surface of a field function f : R 3 → R, defined as:

S = {p ∈ R 3 | f (p) = C}. (1)
We can identify two types of field functions: globally supported and locally supported field functions.

In fields of global support, implicit surfaces are defined as 0-isosurfaces and the inside of the surface can be either defined by points p in space such that f (p) < 0 as in Hoffman et al. [1985] and [START_REF] Barthe | Extrusion of 1d implicit profiles: Theory and first application[END_REF] works, or by points p in space such that f (p) > 0 as in [START_REF] Pasko | Function representation in geometric modeling: concepts, implementation and applications[END_REF] works.

Local support field functions are commonly positive, decreasing functions of the distance, for instance, to a skeleton, that vanish at a distance R from this skeleton. A skeleton can be a point, a line segment, etc. In this case, the implicit surface is an iso-surface "around" the skeleton, e.g. a sphere if the skeleton is a point. The implicit surface is the 0.5-iso-surface, and we use the convention that f (p) > 0.5 defines the interior and f decreases outside an object.

Composing field functions: One of the strengths of implicit modeling is the ability to combine fields to form a new shape, using a composition operator as:

f = g(f 1 , f 2 ) (2)
where g : R 2 → R is a binary composition operator, f 1 and f 2 are field functions defining the combined implicit surfaces and f is the field function defining the resulting object. As we will see, different operators are to be developed for global and local field functions.

Blending: A standard operator for generating smooth blends is the sum [START_REF] Blinn | A generalization of algebraic surface drawing[END_REF]]:

g(f 1 , f 2 ) = f 1 + f 2 . (3)
Unfortunately, it does not provide any control on the shape of the resulting blend, and can only be applied to local support field functions or to convolution surfaces. Therefore, several binary blending operators were proposed over the years for field functions with global support [Bloomenthal 1997b]. The idea is to blend 0-isosurfaces while meeting specific limit properties, continuity constraints or aspect of iso-curves [Hoffmann and Hopcroft 1985;[START_REF] Rockwood | The displacement method for implicit blending surfaces in solid models[END_REF][START_REF] Barthe | Extrusion of 1d implicit profiles: Theory and first application[END_REF]. For instance, the superelliptic blend [START_REF] Rockwood | The displacement method for implicit blending surfaces in solid models[END_REF]] is defined as:

g(f 1 , f 2 ) = 1 - 1 -f 1 r 1 t + - 1 -f 2 r 2 t + (4)
where [ * ] + = max(0, * ), t controls the shape of the blend, making it closer to the combined objects when it increases, r 1 and r 2 control the boundaries of the blend on each combined primitive. This operator exhibits some gradient discontinuities that can be avoided by using Pasko et al. blending set theoretic operator [1995] defined as:

g(f 1 , f 2 ) = f 1 + f 2 -f 2 1 + f 2 2 + a 0 1 + f 1 a 1 2 + f 2 a 2 2 (5)
where a 0 , a 1 and a 2 have an equivalent effect on the blend than, respectively, t, r 1 and r 2 in Equation 4. Better field variations and user control can be achieved using Barthe et al. functional blend [2001] defined as:

g(f 1 , f 2 ) = min(f 1 , f 2 ) -H (|f 1 -f 2 |) (6)
where H : R → R is a piecewise cubic polynomial defining the shape of the blend from user selected control points.

Union: Another well known operator performs a union using the max function [START_REF] Sabin | Function representation of solids reconstructed from scattered surface points and contours[END_REF][START_REF] Ricci | A Constructive Geometry for Computer Graphics[END_REF]:

g(f 1 , f 2 ) = max(f 1 , f 2 ). (7)
This operator generates C 0 continuous field functions that can cause creases on subsequent blends with additional surfaces.

Clean union: Clean union operators have been introduced by [START_REF] Pasko | Function representation in geometric modeling: concepts, implementation and applications[END_REF] (Equation 8) and improved by [START_REF] Barthe | Two-dimensional potential fields for advanced implicit modeling operators[END_REF] for field functions with global support.

g(f 1 , f 2 ) = f 1 + f 2 -f 2 1 + f 2 2 (8)
They improve on the standard max operator whose sharp field is illustrated in Figure 2(a-bottom) by performing a union of the surfaces while providing a smooth field function everywhere else ). In a sense, a clean union operator is a union operator on the combined surfaces and a blending operator of their O. Gourmel et al.

field functions everywhere else. This allows an object built using such an operator to be later blended with another implicit primitive without introducing unwanted sharp edges into the blend [START_REF] Pasko | Function-based shape modeling: mathematical framework and specialized language[END_REF]. plots of g(f 1 , f 2 ) where the axes are f 1 and f 2 . Black lines: some isocurves, white line: 0.5-iso-curve (surface of interest), background color: gradient direction with vertical (resp. horizontal) gradient in blue (resp in red), green lines: boundaries outside of which g = max(f 1 , f 2 ). Second row: implicit surface g(f 1 , f 2 ) = 0.5. Some iso-lines of g are depicted in a horizontal cutting plane.

Visual representation of binary operators:

In their work, Hoffman and Hopcroft [1985] introduce an R 2 space in which operators g are defined and plotted with the value of the field functions f 1 and f 2 respectively as abscissa and ordinate. In this space, vertical lines represent the iso-surfaces of the field function f 1 and the horizontal lines, those of the field function f 2 . The top row of Figure 2 illustrates three operators plotted in this space in which the white line represents the object iso-surface. The composition operator g = max is represented as in Figure 2(a) while a clean union operator is represented as in Figure 2(b) and a blending operator that smoothly links iso-surfaces of the field function f 1 with those of the field function f 2 at the vicinity of their intersection (i.e. close to

f 1 = f 2 ) is represented as in Figure 2(c).
Composing locally supported field functions: When using locally supported field functions, primitives are defined by 0.5-iso-surfaces and field functions uniformly equal zero outside the boundary of their support. Thus, composition operators have to be such that g(0, f 2 ) = f 2 (on and outside the support of f 1 , g reproduces the field function f 2 ) and g(f 1 , 0) = f 1 (on and outside the support of f 2 , g reproduces the field function f 1 ). This also implies that g(0, 0) = 0. These properties make the definition of composition operators more subtle and tedious than for global supports. All these properties are fulfilled by the operators presented in Figure 2. For instance, the clean union operator presented in [START_REF] Barthe | Controllable binary csg operators for "soft objects[END_REF], similar to the one depicted in Figure 2(b), is defined as:

g(f 1 , f 2 ) =              max(f 1 , f 2 ) if f 1 ≤ 1 2 and f 2 ≥ f 1 2 or f 2 ≤ 2f 2 1 or f 1 > 1 2 and f 2 ≤ f 1 2 or f 2 ≥ 2f 2 1 {C : ĥC (f 1 , f 2 ) = 1} otherwise (9)
in which ĥC is defined as:

ĥC (f 1 , f 2 ) =            (f1-2C 2 ) 2 +(f 2 -2C 2 )) 2 C-2C 2 if f 1 ≤ 1 2 f 1 - √ C 2 2 + f 2 - √ C 2 ) 2 C- √ C 2
otherwise.

(10)

An intricate closed form solution for the computation of C is given in [START_REF] Barthe | Controllable binary csg operators for "soft objects[END_REF]]. For blending operators, the sharpness of the blending shape can be controlled by a parameter t if the sum is replaced by [START_REF] Ricci | A Constructive Geometry for Computer Graphics[END_REF]]:

g(f 1 , f 2 ) = f t 1 + f t 2 1 t . (11)
However, if some additional control is required on the blend boundary, as was achieved in the global support case through Equations 4, 5 and 6, then the blending operator proposed in [START_REF] Barthe | Controllable binary csg operators for "soft objects[END_REF] needs to be used. This operator is illustrated in Figure 2(c) with θ = θ 1 = π 2 -θ 2 and defined as:

g(f 1 , f 2 ) =    max(f 1 , f 2 ) if f 2 ≤ tan(θ 1 )f 1 or f 1 ≤ cot(θ 2 )f 2 {C : hC (f 1 , f 2 ) = 1} otherwise (12)
in which f 2 = tan(θ 1 )f 1 and f 1 = cot(θ 2 )f 2 correspond to the green lines in Figure 2(c) (the max is used under the first line and on the left of the second). In between these two lines, hC is defined as:

hC (f 1 , f 2 ) = (f 1 -C. cot(θ 2 )) 2 (C -C. cot(θ 2 )) 2 + (f 2 -C. tan(θ 1 )) 2 (C -C. tan(θ 1 )) 2 (13)
and θ 1 , θ 2 control the boundaries of the blend on each combined primitive (as r 1 , r 2 in Equation 4and a 1 , a 2 in Equation 5).

Opening angles to control the blend: Our solution is based on the work by [START_REF] Barthe | Controllable binary csg operators for "soft objects[END_REF] we already mentioned, in which the operator g, is designed to:

(1) give fine control over the resulting shape, as does the superelliptic blend [Hoffmann and Hopcroft 1985;[START_REF] Hsu | Field functions for blending range controls on soft objects[END_REF],

(2) avoid the gradient discontinuities in the field function exhibited by displacement-blends [START_REF] Rockwood | The displacement method for implicit blending surfaces in solid models[END_REF]],

(3) handle field functions with local support.

The solution presented in Equation 12uses the concept of an opening angle θ ∈ [0, π/4] (figure 2(c)) first introduced in [START_REF] Barthe | Two-dimensional potential fields for advanced implicit modeling operators[END_REF]]. The blend is delimited by boundary lines (green lines in figure 2(c)). Inside these boundaries the iso-lines of g(f 1 , f 2 ) are arcs of ellipses defined in Equation 13(as suggested in [Hoffmann and Hopcroft 1985]). Outside, no blend occurs and the resulting field function is defined as In practice, the opening angle can be set automatically from a point p, selected in object space R 3 , on one of the primitives where the user wants to switch between blending and union. For instance, suppose that the point p(x, y, z) is selected on the 0.5iso-surface of the field function f 1 , i. Interpolating between clean union and blending: Operators interpolating between clean union and blending are able to perform both union and a smooth blend in the same composition. For instance, they enable the combined objects to blend where they are in contact while avoiding any shape deformation everywhere else.

f = g(f 1 , f 2 ) = max(f 1 , f 2 ) (Equa- tion

GRADIENT-BASED COMPOSITION

Let f 1 and f 2 be the field functions of the implicit primitives to be composed and g be the composition operator we are looking for.

As in [START_REF] Rockwood | The displacement method for implicit blending surfaces in solid models[END_REF]], we claim that g should not only depend on f 1 and f 2 , but also on their gradients. Indeed, this appears as a pertinent option to smoothly blend intersecting shapes where they form a sharp angle, while avoiding bulges where their surfaces are already aligned. This section first presents the key features of our solution. We then explain the field continuity issues that motivate the introduction of a family of quasi-C ∞ operators in Section 4.

Controlling opening from angle between gradients

As we just mentioned, g should be able to generate smooth blending between selected parts of the input shapes and union elsewhere.

Following [START_REF] Barthe | Controllable binary csg operators for "soft objects[END_REF]] we choose g with an opening angle parameter θ ∈ [0, π/4] (see Figure 2(c) and Section 2). θ controls the locality of the blend, i.e. the limit values in the (f 1 , f 2 ) space out of which a union is computed.

The key feature of gradient-based blending is to define the opening angle θ as a function of the angle α between the field gradients at the query point p:

θ = θ(α(p)) with α(p) = arccos ∇f 1 (p) • ∇f 2 (p) ∇f 1 (p) ∇f 2 (p) ( 14
)
where α is the angle between the field gradients and the opening angle θ : [0, π] → [0, π/4] is defined as a continuous function we call an opening function.

While a constant opening function θ applies the same type of composition everywhere along the input surfaces (such as smooth blending, shown in the first row of Figure 4), choosing opening functions that adequately tune the opening angle according to α allows us to seamlessly solve the four challenging problems we mentioned. Figure 4 shows how this is achieved; from left to right we have:

(a) The bulging problem: the 'T' in the top row shows the inflation produced by smooth blending where surfaces are aligned (α close to 0). To avoid this, the opening function should switch to union when α = 0 (aligned gradients). This is done by choosing θ such that θ(0) = π/4 (union for aligned gradients) and θ(π/2) = 0 (blending for orthogonal gradients).

(b) This behavior also naturally prevents the absorption problem on small details since the inflation produced by the blend is gradually reduced as we come close to the "top" of sharp features (where the input surfaces are aligned).

(c) The locality problem occurs when two surfaces come close to each other and blend while they do not intersect, so when α is close to π. This is easily avoided by choosing θ such as θ(π) = π/4 (union for opposing gradients) and θ(π/2) = 0 (blending for orthogonal gradients).

(d) This setting also avoids the topology problem, i.e. the change of topological genus shown on top, since there are points inside each handle where the angle between gradients is π and which will therefore remain outside of the composed shape. The opening functions used to generate these compositions are depicted in Figure 4(e). Note that the one on the second row, which we call the Camel opening function for its shape, solves all four challenging problems we mentioned. It will be fully described in Section 5. The constraint for g and θ to be C ∞ is in fact theoretical: in practice, we just need g and θ to approximate some C ∞ functions at a given precision. This has to be done without introducing oscillations in the field, since the latter would result in curvature artifacts. We say that a function verifying this looser constraint is quasi-C ∞ continuous. Relaxing our continuity constraint to quasi-C ∞ will enable us to efficiently store our operators as textures on the GPU and to evaluate them with hardware linear interpolation.

Continuity issues

Section 4 introduces our opening functions and our new composition operators, which are applied to shape modeling and animation in Section 5. In these sections, the input fields are supposed to be singularity free, i.e. with smoothly varying, non-zero gradients vectors. The application of our framework to more general settings will be discussed in Section 6.

QUASI-C ∞ OPERATORS 4.1 Opening functions

In this section we describe a practical set of C ∞ functions for defining the opening functions θ that convert the angle between gradients into the opening angle of the composition operator. Based on the work of section 3 this family of functions should be flexible enough to implement various modes of conversion, including, but not restricted to those of Figure 4(e). Opening functions designed to model organic shapes and to animate contact situations will be shown in Section 5. The family of functions needed for θ should enable to set some key values and provide some slope control; without which, unwanted bumps may appear, depending on the slopes of the input fields around regions where a clean union is applied. See Figure 6. We therefore use a class of functions θ(α), α ∈ [0, π] defined from eight parameters α 0 , α 1 , α 2 , θ 0 , θ 1 , θ 2 , w 0 and w 1 :

θ(α) =          θ 0 α ≤ α 0 κ α-α 1 α 0 -α 1 w 0 (θ 0 -θ 1 ) + θ 1 if α ∈]α 0 , α 1 ] κ α-α 1 α 2 -α 1 w 1 (θ 2 -θ 1 ) + θ 1 if α ∈]α 1 , α 2 [ θ 2 otherwise ( 15
)
where

κ(x) = 1 -exp 1 - 1 1 -exp 1 -1 x . ( 16
)
This function is controlled with eight parameters: three values of α i (i = 0, 1, 2), where α 0 and α 2 are limits of the intervals [0, α 0 ] and [α 2 , π] in which θ is set to be constant, and the three associated values for θ (θ 0 = θ(α 0 ), θ 1 = θ(α 1 ), θ 2 = θ(α 2 )) and two additional parameters w 0 and w 1 controlling the slope of the opening function in respectively [α 0 , α 1 ] and [α 1 , α 2 ]. See Figure 7. In order to ensure C ∞ continuity, θ is piecewise defined from a function κ of class C ∞ that is flat at the boundaries of its support (i.e. all their successive derivatives equal zero).

θ 1 θ 2 θ 0 α 0 α 1 α 2 π 0 π/4 α 0 =0 α 1 =π/2 α 2 =π θ (a) (b)

Composition operators

Let us now define another family of quasi-C ∞ functions for the composition operator g. This operator should produce smooth blends and also be able to generate other useful compositions, such as bulges when soft objects come in contact [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF]]. Whatever the type of composition, the operator should be parameterized by the opening angle θ, θ ∈ [0, π/4]. The challenge is to define an operator that will fully apply the composition when θ = 0 (opened operator, e.g. a smooth blend), and gradually switch to a clean union between the input shapes when θ = π 4 (closed operator), while avoiding the limitations of Bernard's operator (Section 1.1 and Figure 3).

Following the construction of operators with an opening angle, our operator g uses the angle θ to define the width of the region where the composition is applied at the level of the iso-value (red isocurve in Figure 8). We also use boundary curves k θ bounding the region where the values of g modifies the field of the combined field functions f 1 and f 2 (light grey area in Figure 8). In this region, g = ḡ and g = max outside. A critical contribution here is the special shape of our boundary curves k θ when θ varies from 0 to π/4 (Figure 8 These boundary curves are defined by C ∞ functions k θ (f ), where f is either f 1 or f 2 . We set k θ (f ) = tan(θ)/2 inside of the implicit primitives, where f ≥ 0.5. Outside the primitives, where f < 0.5 the expression of k θ is:

k θ (f ) = tan(θ) 2 4 1 + tan(θ) λ θ (f ) 2 , ( 17
)
with and φ is a function ensuring the C ∞ continuity of λ θ . φ is computed by binary search at the required precision, from φ -1 defined as:

λ θ (f ) = f if f ≤ tan(θ) 2 1-tan(θ) 4 φ(2 2f -tan(θ) 1-tan(θ) ) + tan(θ) 2 otherwise ( 18
) • O. Gourmel et al. φ s(φ) 1 1 k θ f 1 f 2 φ (f , f ) 1 2 C k θ C ( ) C k θ C ( ) C -kθ C ( ) C -kθ C ( ) ) ( s ( ) φ - O k θ f 1 f 2 (a) (b) (c)
φ -1 (x) = x + 1 e log log 1 ν(x) + 1 + 1 , ( 19
)
where e = exp(1) and with

ν(x) = exp (exp(e -ex) -1) -1. (20)
In practice, the computation of φ does not affect efficiency since our operators are pre-computed on the GPU (see Section 4.3). The value of g inside the region bounded by the boundary curves k θ (in grey on Figures 8 and9 (c)) is defined by an operator ḡ : R 2 → R, which sets the type of composition we desire. The equation of our composition operator g is thus:

g(f 1 , f 2 ) = max(f 1 , f 2 ) if f 1 ≤ k θ (f 2 ) or f 2 ≤ k θ (f 1 ) ḡ(f 1 , f 2 ) otherwise (21)
The function ḡ is defined from a silhouette curve s(ϕ) : [0, π/2] → R describing the shape of the iso-curves of ḡ as illustrated in Figure 9. For instance, an arc-of-a-circle is created using s(ϕ) = 1, ∀ϕ. For given values of f 1 and f 2 , the value of our operator ḡ(f 1 , f 2 ) is computed by solving the following equation in C:

ḡ(f 1 , f 2 ) = {C : hC (f 1 , f 2 ) = 1} (22) with hC (f 1 , f 2 ) = (f 1 -k θ (C)) 2 + (f 2 -k θ (C)) 2 s(ϕ) (C -k θ (C)) (23) and ϕ = arctan f 2 -k θ (C) f 1 -k θ (C) (24)
which expresses the fact that the iso-curves of ḡ are scaled versions of s. In practice, equation ( 23) does not always have a closed form solution. We thus precompute sampled values of g and interpolate them on the GPU, as detailed in Section 4.3.

Using any C ∞ silhouette curve s, such that s(0) = s(π/4) = 1 and with flat derivatives at the level of the boundary curve (so that g's derivatives match those of the max function) is thus sufficient for getting a C ∞ operator g. To illustrate how different these silhouette curves can be and the variety of resulting compositions, Section 5 will give the equations of the curves we use for smooth blending and for a "bulge-in-contact" composition.

Implementation on the GPU

As discussed above, we use families of opening functions and of composition operators that have all the required smoothness properties, but that may not have any closed form expression. For efficiency, we developed a GPU implementation for these functions, which enables us to apply gradient-based compositions in realtime. The GPU implementation is based on the following precomputations: the opening function which gives θ from the angle between gradients (equation 15) and its derivative ∂θ ∂α are stored as 1D textures. We also pre-compute values of the operator g and of its partial derivatives ∂g ∂f 1 , ∂g ∂f 2 , ∂g ∂θ as functions of θ, f 1 and f 2 in a regular grid stored as a 3D texture. This is done efficiently using the following procedure: The values of g are computed by grid slice (f i 1 , f j 2 ), (i, j) ∈ [0; N -1] 2 , for each grid coordinate θ k , k ∈ [0; N -1] (N 3 being the size of the grid). In a slice, for each i, we set C i = f i 1 and we follow the C i -iso-curve G i of g and update the sample values as follows (see Figure 10): -Outside of the boundary curves, i.e. for each j such that f j 2 < k θ (f i 1 ), g = max(f 1 , f 2 ) (Equation 21), so we set g(f i 1 , f j 2 ) = C i and g(f j 2 , f i 1 ) = C i . -Inside of the boundary curves, for each sampling point X(X x , X y ) located in the bounding box of the silhouette curve s scaled by (C i+1 -k θ (C i )) with the bottom-left corner in O i , we compute γ 1 and γ 2 as in Equations ( 26) and ( 27). When γ 1 ≥ 0 and γ 2 ≥ 0, X is between G i and G i+1 , i.e. the C i -iso-curve and the C i+1 -iso-curve, and we compute ḡ(X) by linearly interpolating C i and C i+1 as in Equation ( 25).

X k θ φ i O i G i+1 k ( θ C ) i C - i γ 1 G i f 1 f 2 C i C i+1 X k θ φ i+1 O i+1 G i+1 γ 2 G i f 1 f 2 C i C i+1 k ( θ C ) i+1 C - i+1 (
ḡ(X) = γ 2 C i + γ 1 C i+1 γ 1 + γ 2 (25)
with

γ 1 = X -O i -s(ϕ i ) (C i -k θ (C i )) (26) γ 2 = s(ϕ i+1 ) (C i+1 -k θ (C i+1 )) -X -O i+1 (27)
where

O i = (k θ (C i ), k θ (C i )) and ϕ i = arctan X y -k θ (C i ) X x -k θ (C i ) . (28)
Once all the values of g are computed, its partial derivatives are evaluated using finite differences and stored in other grids.

During field queries for a composed implicit surface, the values in these textures are linearly interpolated on the GPU: we first get θ from the gradients of the input fields and then g(f 1 , f 2 ) with this specific opening angle θ. Depending on the primitive, gradients are either computed using a closed form expression or finite differences. To efficiently apply successive compositions, the gradient of g(f 1 , f 2 ) is computed as:

∇g = ∇f 1 ∂g ∂f 1 + ∇f 2 ∂g ∂f 2 + ∇α ∂θ ∂α ∂g ∂θ (29)
A consequence of our GPU implementation is that the operators we use in practice are not exactly the C ∞ functions we developed, but quasi-C ∞ approximations of the latter: indeed, linear interpolation of values in a texture prevents unwanted oscillation and enables us to approximate the functions at any desired precision. We tested the method with 16 3 , 32 3 , 64 3 and 128 3 grids for g and its gradients. Figure 11 shows the improvement of reflection lines on a very close view of a composed primitive when texture resolution increases. Our current implementation uses 128 3 grids which is accurate enough to prevent any visual artifact, even after multiple compositions. Note that if a very high precision was required, e.g. for a manufacturing application, an off-line evaluation mechanism with bounded error could be set up, where binary search would be used to find the resolution needed to obtain a given error bound.

APPLICATIONS

In this section, we present applications of our generic composition operators to constructive modeling and animation. This leads us to define specific silhouette curves s and the associated opening functions θ.

16 3 32 3 64 3 128 3 max error 5.03 × 10 -2 2.90 × 10 -2 1.39 × 10 -2 6.72 × 10 -3 mean error 2.97 × 10 -3 1.43 × 10 -3 7.32 × 10 -4 3.69 × 10 -4

Fig. 11. Close-up on the surface obtained after three successive blends when g is stored into textures of different resolutions. The piecewise linear interpolation generates artifacts in reflections (here illustrating the continuity of the fifth derivatives) with the 16 3 and 32 3 grids that are not perceptible using higher resolutions. The third and fourth rows respectively give the max and mean approximation error provided by the tri-linear interpolation of the pre-computed values.

Smooth blending for constructive modeling

Being able to smoothly blend arbitrary shapes is one of the most useful compositions in constructive modeling systems. To create a good blending operator ḡb , we need a C ∞ silhouette curve sb , which, in addition to meeting the value and derivative constraints listed at the end of Section 4.2 (values 0.5 and null derivatives of every order at the extremities), avoids curvature oscillations as much as possible, as illustrated in Figure 12.

Our solution plotted in Figure 12(c) is a symmetric silhouette curve whose curvature monotonically increases, with a maximum at ϕ = π/4, then decreases back to zero. We define sb as follows:

sb (ϕ) = 1 - 1 e log log 1 ν(ϕ) + 1 + 1 (30) • O. Gourmel et al.
Barthe's blending "camel" blending "organic" blending "organic" blending in interactive modeling Fig. 14. Comparison of Barthe's blending with our operator g b , used either with the "camel" or with the "organic" opening function.

where ν is defined in Equation 20.

For our pre-computations of g on the GPU, we express this profile function in polar form, as sb (ϕ), by sampling the curve at points X i = (x i , sb (x i )) and pre-computing the value of ϕ i at each X i .

To evaluate the function sb (ϕ), we use binary search to find the interval [ϕ i ; ϕ i+1 ] containing ϕ, and then use linear interpolation.

Although not exactly flat at the boundaries, this silhouette curve has zero derivatives in ϕ = 0 and ϕ = π/2 up to an error of ǫ = 10 -5 , which is accurate enough for our needs. This operator is illustrated in R 2 space in Figure 13 This blending operator g b computed using sb can be combined with different opening functions according to the need. Our applications to shape modeling lead us to design two specific examples of opening functions.

The first one, whose parameter values are given in the first row of Table I, is called the camel opening function because of the shape of the associated curve (see its plot in Figure 15(a)). It is symmetric and was especially designed to solve, for Wyvill's local-support soft-objects [START_REF] Wyvill | Data Structure for Soft Objects[END_REF]], all four problems mentioned in the introduction (Figure 4). In the case of CAGD applications, where two primitives join, this opening function provides smooth, predictable blending behavior while preventing bulges.

The second one, whose parameter values are given in the second row of Table I, is called the organic opening function (see its plot in Figure 15(b)). It was developed for an interactive modeling system dedicated to the modeling of organic shapes, using globalsupport convolution surfaces based on Cauchy kernel [START_REF] Tai | Prototype modeling from sketched silhouettes based on convolution surfaces[END_REF]]. Here, being able to model both smooth and sharp features in the same composition is important: for instance, when blending a character's nose to the head (Figure 1(b)) or the leg and tail of an animal to its body (Figure 14 (right)), the junction has to be smooth at the top but sharp at the bottom. To achieve this, our organic opening function is set to be asymmetric: surfaces smoothly blend where they have nearly orthogonal gradients, while a clean union generating sharp creases is used when gradients tend to be in opposite directions (see Figure 14(left)). This enables us to get detailed organic-like models with both smooth parts and sharp creases.

Table I. Parameter values defining our different opening functions.

Bulge in contact for soft objects animation

To illustrate the variety of compositions we can achieve, we designed a "bulge-in-contact" operator g c , inspired by [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF]], which instead of blending primitives that intersect, mimics the effect of two surfaces in contact that bulge as if made from a soft material such as rubber (Figures 16, 1 Keeping in mind that the silhouette curve s directly gives the shape of the deformation, Figure 16(a) shows the deformation of f 1 on the right and of f 2 on top. We set s such as to create smoothly increasing and then decreasing bulges. This is done by using values smaller than 0.5 on the silhouette curve. Defining these contact silhouette curves is done using the parametric expression sc (t)(x c (t), y c (t)) with:

x c (t) = t if t ≥ 1 1 -K exp 1 -1 2-t (t -1) otherwise y c (t) = 1 -K exp 1 -1 t (1 -t) if t ≥ 1 2 -t otherwise (31)
for t ∈ [0, 2], where K ∈ [0, 1] is a parameter that controls the thickness of the bump. The effect of K on the iso-curves of ḡc is illustrated in Figure 16 (a). The curve sc is perfectly flat at the extremities and is, as desired, only C 0 continuous at ϕ = π 4 in order to simulate the fold to be created between the two input surfaces. Again, the polar form sc (θ) can easily be evaluated using binary search, as done for sb in Section 5.1.

Finally, a contact opening function, whose parameter values are given in the third row of Table I, is designed so that when used with our bulge-in-contact operator, objects do not deform at a distance but bulge when they intersect, mimicking contact between soft objects, as illustrated in Figure 15(c). The operator g c is illustrated with K = 0.8 in Figures 1(d) and 16(c).

RESULTS AND DISCUSSION

Even though the implementation of the set of functions presented here can be considered complex, its usage is not. The functions can be used in an interactive modeling framework as any other operator; objects are added to the scene, selected and the operator is applied to two of them. With most traditional implicit operators the resulting shape is fixed after applying the operator. In contrast our new set of functions provide ways to further modify the shape, by altering the silhouette curve and the opening function in a graphical editor, however, most of the useful opening functions are defined in this paper and cover most of the practical cases. New opening functions only need be defined for specialized operations.

The usefulness of our generic composition framework for constructive modeling and animation applications is illustrated in Figures 17 and 1(d). The model presented in Figure 17 has been built using ten different primitives combined with either our blending operator together with the camel opening function or with our clean union operator (generating sharp edges). In the case of a blend overlapping a sharp edge, we use the procedure presented in Section 7 and in Figure 19. The object has been created interactively using the Blob-tree [START_REF] Wyvill | Extending the csg treewarping, blending and boolean operations in an implicit surface modeling system[END_REF] in several sessions over a period of about eight hours. In the "lava lamp" animation, g c is used in conjunction with our blending operator to also get some blending when bubbles leave the bottom of the lamp. Other interesting animation effects can be easily produced by dynamically modifying the opening function parameters during an animation. For instance, Figure 18 illustrates an animation where primitives do not blend when they are not in contact, but do blend (and thus deform each other) when they separate even though they are not in contact anymore, as in standard blending. Achieving good quality results leads us to use more complex equations for composition than in all previous methods. The efficient evaluation of our operators is based on pre-computations that only need to be performed once for each operator and for each opening function, and just requires the memory space of a few textures. This implementation allows us to use arbitrarily complex expressions for g and θ, while answering field queries at the cost of a simple texture fetch. Therefore, when designing a new operator, one can focus on the desired effect rather than on the best properties a closed form equation could achieve. The precomputation of our operator and its partial derivatives in 128 3 grids takes 600 ms on a Core I7 950 and their transfer to the device memory from the host memory takes 3 ms. The 786 million evaluations of g required for the rendering of Figure 11 take less than 270 ms on a NVIDIA GTX 480.

LIMITATIONS

Firstly, as we emphasized in Section 4.1, adequate slope for the opening function has to be used, depending on the slope of the input field functions, for not creating extra artifacts when removing a bulge. In practice, we use the interactive feedback of implicit modeling systems to pre-set the steepness parameters w 0 and w 1 of our opening functions according to the family of field function used. Typically, we need sharper opening functions to blend softobjects based on smooth polynomial field functions than to compose sharper convolution surfaces. Once the steepness parameters are set, we design specific opening functions by adjusting the three input values for (α i , θ(α i )). Providing an automatic tuning of the opening function's slope would be a better solution. The method, as currently implemented, still requires little user interaction since the manual tuning only needs to be done once for a given family of field functions, for each opening function.

A drawback of our gradient-based operators is that they theoretically require at least the continuity of the gradient directions of the field functions f 1 and f 2 to be combined. Using our composition framework in a more general setting, where the field may have a few singular points (with null gradients and therefore undefined angle between gradients) would however be desirable, as well as using it even for shapes with gradient discontinuities, e.g. for primitives with sharp edges. In case of null gradient vectors, we simply decrease the opening angle to min(θ 0 , θ 1 , θ 2 ) when the gradient norm becomes close to null, so that the opening function is constant at the vicinity of the singularity. This only leaves us with the case of artifacts created by discontinuities of the gradient field. This is illustrated by Figure 19(left), where the camel opening function is used for blending in the vicinity of a sharp edge. The problem can be solved in the following way: while primitives are not in contact, an opening function preventing blend is required, and as soon as the primitives intersect, a blend can be progressively performed. This is done easily by manually tuning a few opening function parameters as illustrated in Figure 19. In an animation, the opening function parameters have thus to be pre-set in key primitive positions such as those presented in Figures 19 andthen As with all recent composition methods, the gradient-based operators are restricted to binary compositions of implicit surfaces. This may limit their applicability in situations where many primitives need to be symmetrically connected at once. Generalizing the idea of using gradients into n-ary composition operators would be a topic for future work. Similarly, although they can be used to compose convolution primitives, our operators cannot replace the unwanted handle (camel) clean union while primitives are not in contact blending when primitives intersect Fig. 19. Composition of a cube and a sphere at a gradient discontinuity (an edge of the cube). On the left, the artifact generated when the camel opening function is used for blending, and then, opening functions solving the problem when the sphere moves and intersects the cube.

sum performing the convolution integral of each element (points, segment, triangles) of a complex skeleton defining a convolution surface [Bloomenthal 1997b]. We do not provide a solution to the problem of a single convolution primitive, from unintentionally blending with itself.

CONCLUSION AND FUTURE WORK

We have presented a generic family of composition operators, applicable to both constructive implicit modeling and animation. These operators, implemented in real-time on the GPU, allow precise tuning of the transition between smooth blends and sharp union within a single composition operation, making compositions both more general and more predictable. The main feature of our approach is the use of gradient-controlled blending: our operators are not only a function of the field functions to be combined, but also of their gradient. This avoids the well known weaknesses of implicit modeling, such as creating unwanted bulges and blending at distance. In general, the topology of the resulting shape can be set to the topology of the union. Moreover, our operators naturally prevent small details from blowing up even when blended into much larger primitives, thanks to the sharp changes of gradient values.

In consequence, complex models with both smooth parts and sharp creases can be easily created, in contrast with the common idea that implicit surfaces can only represent smooth, blobby shapes.

The natural extension of using the angle between field gradients would be the use of the information on the norm of these gradients to improve blends. This norm could be used to automatically set the slope parameters of the opening functions, according to steepness of the input fields. Another extension would be to rely on the second derivatives of the input fields to detect fast local gradient variations. This would help us to avoid user interaction currently necessary to adjust opening functions in areas where the field functions are very distorted or at the vicinity of gradient discontinuities such as sharp edges. Lastly, the general methodology we developed for the opening functions, i.e. defining C ∞ curves with some shape control parameters, could be re-used to introduce local-support, C ∞ fields for implicit primitives. This would ensure always having smooth gradient fields, whatever the number of gradient-based compositions.

  Fig. 1. (a) The unwanted bulge created by standard blending (left) is suppressed by our gradient-based method while preserving the desired topological genus (right). (b) We also prevent small details from blowing up and primitives from blending where they do not intersect. (c) Our generic framework allows the modeling of a wide range of composition operators such as "bulge in contact" and (d) blending (left) and "bulge in contact" (right) combined in the same animation to mimic a lavalamp.

Fig. 2 .

 2 Fig.2. Three compositions applied to a pair of cylindrical implicit primitives forming a cross: (a) standard union based on max, (b) clean union, (c) Barthe's blending operator parameterized by an opening angle θ. First row: plots of g(f 1 , f 2 ) where the axes are f 1 and f 2 . Black lines: some isocurves, white line: 0.5-iso-curve (surface of interest), background color: gradient direction with vertical (resp. horizontal) gradient in blue (resp in red), green lines: boundaries outside of which g = max(f 1 , f 2 ). Second row: implicit surface g(f 1 , f 2 ) = 0.5. Some iso-lines of g are depicted in a horizontal cutting plane.

  12) in order to exactly reproduce one of the fields of the input field functions. The boundary lines are controlled by varying the opening angle θ so that any intermediate blending (illustrated with θ = π/8 in figure2(c)) can be produced between a full blend when θ = 0 and a union (g = max) when θ = π/4 (figure2(a)). When θ ∈ [0, π/4], this operator is G 1 continuous and when θ = π/4, it is C 0 .

  Fig. 3. Bernhardt's operator interpolating between clean union and Barthe's blending. The distortions of the iso-curves of g (a) result in a G 1only implicit surface (b) which still has a bulge on top and exhibits a poor curvature distribution, shown by depicting reflection lines (c).

  They have been introduced by[START_REF] Pasko | Bounded blending for Function-Based shape modeling[END_REF] for globally supported field functions. Pasko's operator is defined by Equation 5 in which the coefficient a 0 varies: if a 0 = 0, the operator performs a clean union (as in equation 8), and the larger a 0 , the larger the blending. Bernardt at al.[2010] proposed a such operator for locally supported field functions. This operator linearly interpolates between the clean union operator defined in Equation 9 and the blending operator defined in Equation12. The limitations of these methods have been discussed in Section 1.1 and Bernhardt's operator is illustrated in Figure3. Note that the clean union operator introduced in[START_REF] Barthe | Controllable binary csg operators for "soft objects[END_REF]] (Figure2(b)) does not support the insertion of an opening angle that would allow the interpolation between a blend and a clean union.

•O.

  Gourmel et al.

Fig. 4 .

 4 Fig. 4. Comparison between standard blending (first row) -equivalent to a null opening function -and our method (second row) with a opening function that smoothly switches from clean union when α = 0 or α = π to blending when α = π/2. Our solutions to: (a) unwanted bulge, (b) blow up of small details, (c) blending at distance, (d) topology modification. (e) Plot of θ(α).

  Fig. 5. (a) Two cylindrical implicit primitives are combined using Barthe's G 1 operator, where our gradient-based opening function is used to set the opening angle. The resulting shape still looks smooth. (b) A third primitive is added with the same blending operator. A C 0 sharp crease now appears. (c) Same scene using our new quasi-C ∞ blending operator.Computing composition based on the gradients of the input fields having a finite level of continuity locally causes a loss of one order of continuity for the resulting field. This decrease cannot be avoided, but efforts can be made such that neither the opening function θ nor the composition operator g (based on it) introduce any extra source of discontinuity, which could cause artifacts in case of consecutive blends. In particular, Barthe's operator[START_REF] Barthe | Controllable binary csg operators for "soft objects[END_REF], depicted in Figure2(c), cannot be used in our framework since it is G 1 only and sets up a transition to a standard, C 0 only union. If the opening angle of this operator is set to a function θ of the field gradients, it generates a G 1 field function with the expected blend after the first composition (Figure5(a)) and the loss of continuity after each composition results in a C 0 field function with a sharp edge in the second blend (Figure 5(b)) while a smooth blend as the one illustrated in (Figure 5(c)) was expected. Higher

Fig. 6 .

 6 Fig. 6. Role of the slope of θ in blending. The red (resp. blue) vectors are gradients of the horizontal (resp. vertical) primitive. (a) Union. (b) Blending with a sharp opening function (w 0 = w 1 = 3): unwanted bumps appear where the surfaces are almost aligned, since θ falls very quickly to 0 there. (c) Tuning the slope of θ according to the type of input primitives while keeping its other parameters unchanged solves the problem. Here, we set w 0 = w 1 = 1, a value pre-selected for soft objects.

Fig. 7 .

 7 Fig. 7. (a) The 8 parameters of a opening function θ with w 0 = w 1 respectively equal to 3, 1 and 1/3 for the red, blue and green curves. (b) Three possible opening functions for the blend in Figure 6, using the same values for w 0 and w 1 .

Fig. 8 .

 8 Fig. 8. Our operator g. First row: The opening angle θ defines the width of the region where composition (here, smooth blend) is applied, at the level of the iso-value (red iso-curve). Second row: The boundary curves k θ bound the region where the values of g are computed with ḡ (light grey area). A max operator is used outside. (a) Opened operator (θ = 0). (b) Intermediate situation. (c) Closed operator (θ = π/4).We introduce a family of C ∞ boundary curves whose intermediate shape is shown in the second row of Figure8(b) and in Figure9 (c). These boundary curves are defined by C ∞ functions k θ (f ), where f is either f 1 or f 2 . We set k θ (f ) = tan(θ)/2 inside of the implicit primitives, where f ≥ 0.5. Outside the primitives, where f < 0.5 the expression of k θ is:

Fig. 9 .

 9 Fig. 9. (a) A silhouette curve s(ϕ), (b) defines the shape of the C-iso-curves of ḡ to produce (c) the corresponding operator g.

Fig. 10 .

 10 Fig.10. Computation of the interpolation parameter γ 1 (a) (respectively γ 2 (b)), defined as the distances from X to G i (respectively from X to G i+1 ), used to compute ḡ at a sample point X in Equations 25, 26 and 27.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Close view of the top of the cross formed by two blended cylindrical primitives: (a) using Bernhardt's blending; (b) using our operator g b with a opening function preventing the unwanted bulge; (c) Plot of sb from which g b is derived.

  for different values of θ. Its fairness leads to a clear improvement compared with Bernhart's blending solution as shown with reflection lines in Figures 12(a) and 12(b).

Fig. 15 .

 15 Fig. 15. Three useful opening functions: (a) camel and (b) organic opening functions used with the blending operator g b ; (c) contact opening function used with the bulge-in-contact operator g c .

  (c) and (d)).

  Fig. 16. (a) Plots of sc for K = 0.25 (red), K = 0.5 (green) and K = 0.75 (blue). (b) The corresponding operator g c when θ = π 8 and K = 0.8. (c) The application of our operator g c .

Fig. 17 .

 17 Fig. 17. Constructive modeling examples built using our new operators. In this model, clean union and clean difference operators have been used in conjunction with gradient-based blending operators. Both sharp edges and subsequent gradient-controlled blends have thus been efficiently performed.

  interpolated during the animation.

•O.

  Gourmel et al.

Fig. 18 .

 18 Fig. 18. From left to right, the opening function is modified to create different blending behaviors during an animation.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

ACKNOWLEDGMENTS

This work has been partially funded by the IM&M project (ANR-11-JS02-007), the Natural Science and Engineering Research Council, Canada, and by the advanced grant EXPRESSIVE from the European Research council.


Zhen Zhou
email: zhouzhenqs@mail.nwpu.edu.cn

Feng Chu
email: feng.chu@iup.univ-evry.frache@nwpu.edu.cn

Ada Che

Saïd Mammar
email: said.mammar@iup.univ-evry.fr

A multi-objective model for the hazardous materials transportation problem based on lane reservation

Keywords: hazardous materials transportation, multi-objective optimization, lane reservation, risk, ε-constraint I

This paper presents an application of the lane reservation strategy in the hazardous materials transportation. Once an accident of hazardous materials transportation happens, its effect is significant. Lane reservation can reduce the hazardous materials transportation risk enormously; however, it will also impact on the normal traffic. The proposed problem is to choose lanes to be reserved on the network and select the path for each hazardous materials shipment among the reserved lanes in order to make a trade-off between the impact on the normal traffic and the total risk of the hazardous materials transportation. We first develop a multi-objective integer linear programming model for the problem with multiple types of hazardous materials and multiple origin-destination pairs. Then we adopt the ε-constraint method to solve the multi-objective model. The performance of the method is evaluated on a large set of randomly generated instances. Computational results show that the ε-constraint method can not only solve the proposed problems within reasonable time but also generate candidate solutions that can shape an efficient Pareto front for the multiobjective optimization problem.

INTRODUCTION

In the countries and areas with the developed industry, large quantities of hazardous substances, such as raw materials, intermediate or final products and wastes, are shipped through the transportation networks by different vehicles, like road, rail and inland waterways. For example, in the United States, approximately 800,000 shipments of hazardous materials are transported each day in bulk and in smaller shipment configurations [START_REF]The role of hazardous material placards in transportation safety and security[END_REF].

Hazardous materials include explosives, gases, flammable liquids and solids, oxidizing substances, poisonous and infectious substances, corrosive substances, and hazardous wastes [START_REF] Verter | A GIS-Based framework for hazardous materials transport risk assessment. Risk analysis[END_REF]. Although rare, accidental releases of hazardous materials do occur during transportation, and these events often have very undesirable consequences, including fatalities [START_REF] Verter | A GIS-Based framework for hazardous materials transport risk assessment. Risk analysis[END_REF]. Many researchers has indicated that the risks related to hazardous materials transportation can be of the same magnitude as those arising from fixed installations according to the accident history [START_REF] Vilchez | Historical analysis, of accidents in chemical plants and in the transportation of hazardous materials[END_REF] [START_REF] Leonelli | New detailed numerical procedures for calculating risk measures in hazardous materials transportation[END_REF]. Sometimes the transportation network goes through areas of high population density and areas which are only crossed by transportation routes without having plants, so that, in case of an accidental spill, a large number of persons could be affected. According to U.S. Department of Transportation statistics, 156,442 hazardous materials transportation accidents occurred from 1995 to 2004, resulting in a total of 221 fatalities and 3143 injuries [5]. In European countries, a good many events occurred during the transportation on roads and railways. Specially, the chemical accident in Seveso, Italy in 1976 caused the contamination of a large population by 2, 3, 7, 8-tetrachlorodibenzo-p-dioxin (TCDD). After twenty years later, possible long-term effects examined through mortality and cancer incidence studies indicates that the effect of hazardous material transportation accident is enormous and far-reaching [START_REF] Pesatori | Cancer incidence in the population exposed to dioxin after the "Seveso accident": twenty years of follow-up[END_REF]. Furthermore, hazardous materials releases during transportation may occur in environmentally sensitive areas exposed to high risk value, such as the natural and historical interesting places.

Historical evidences have shown that the hazardous materials transportation risk is often tremendous and thus need to be taken into account in order to keep it under control and reduce it. Generally, we divide the operations research problems closely related to hazardous materials transportation into two main subproblems, as follows.

1) Hazardous materials transportation network design problem

Kara and Verter [7] [8] gave the definition of hazardous materials transportation network design problem; that is, given an existing road network, the hazardous network design problem involves selecting the road segments that should be closed to hazardous materials transportation so as to minimize total risk. They also addressed a bi-level integer linear programming model for the problem where ones can select minimum risk routs as a subset of the transportation network and the carriers select minimum cost on the available subset. Since then, the hazardous materials transportation network design problem has received attention of other researchers, see for Bianco et al [START_REF] Bianco | A bilevel flow model for hazmat transportation network design[END_REF], Erkut and Gzara [START_REF] Erkut | A bi-level programming application to hazardous material transportation network design[END_REF], Erkut and Alp [START_REF] Erkut | Designing a road network for dangerous goods shipments[END_REF].

2) Location and routing problem for Hazardous materials transportation

In recent three decades, many researchers have focused on the location and routing problem for hazardous materials transportation. As we known, the location and routing problem (LRP) is NP-hard. Problems integrating location and routing decisions pertaining to hazardous materials transportation seem much more complicated, because besides the transportation cost and facility cost in the traditional LRP, at least two types of risk are necessarily to be considered: transport risk and facility risk and models for hazardous materials LRP should be multi-objective by nature [START_REF] Erkut | Chapter 9, Hazardous Materials Transportation[END_REF].

The earliest hazardous materials LRP work was done by Shobrys [START_REF] Shobrys | A model for the selection of shipping routes and storage locations for a hazardous substance[END_REF] on locating the storage facilities and selecting routes for the spent fuel shipments so as to simultaneously minimize the total cost and total transportation risk. Jacobs and Warmerdam [START_REF] Jacobs | Simultaneous routing and siting for hazardous-waste operations[END_REF] modeled the hazardous materials LRP as a continuous network flow problem with minimizing the linear combination of cost and risk. It is likely that List and Mirchandani [START_REF] List | An integrated network planar multiobjective model for routing and siting for hazardous materials and wastes[END_REF] firstly considered risk equity as one of the objectives of their model. Current and Ratick [START_REF] Current | A model to assess risk, equity and efficiency in facility location and transportation of hazardous materials[END_REF] followed them. Stowers and Palekar [START_REF] Stowers | Location models with routing considerations for a single obnoxious facility[END_REF] introduced a bi-objective LRP model with a single facility and a single commodity, whose two objectives were different from others, namely minimizing the total exposure and minimizing the maximum exposure.

All the hazardous materials transportation problems listed above doesn't involve in lane reservation strategy. In this study, we address the problem to select an optimal route for hazardous materials transportation shipments based on lane reservation strategy. The essence of the lane reservation problem is to optimally select the lanes and/or the time intervals to reserve for only special tasks so as to satisfy the condition of time or safety and simultaneously minimize the impact on the normal traffic. Lane reservation strategy on the transportation network is a flexible and economic option for special event or situation, such as great sports meeting and great emergencies in the city contingency management. In fact, lane reservation strategy has been successfully applied in true-life by some large sportive games organizers. For example, John Black [START_REF] Black | Strategic Transport Planning, Demand Analysis of Transport Infrastructure and Transport Services for the 27th Summer, Olympiad Held in Sydney, Australia[END_REF] gave a brief description of strategic planning, master planning and operational planning for the 2000 Sydney Olympics, and mentioned the lane reservation for buses only; Zagorianakos [START_REF] Zagorianakos | Athens 2004 Olympic Games' Transportation Plan: a missed opportunity for Strategic Environmental Assessment (SEA) integration[END_REF] stressed on analyzing the importance and advantage of lane reservation strategy for the 2004 Athens Olympics; In France, one of the lanes of A1 between Charles De Gaulle Airport and St Denis is recently dedicated to taxis and buses only between 7am-10am. It is worthwhile to point out that the studies above put emphasis on lane reservation only at the strategic level. At the tactical level, Yingfeng Wu et al. [START_REF] Wu | Heuristic for lane reservation problem in time constrained transportation[END_REF] firstly considered the lane reservation problem in time constrained transportation network. They presented the mathematical model for the problem, and proposed a heuristic to obtain near optimal solutions. Yunfei Fang et al. [START_REF] Fang | Iterative Algorithm for Lane Reservation Problem on Transportation Network[END_REF] proposed an exact algorithm based on cut-and-solve for the lane reservation problem with multiple origin-destination pairs. Then they [START_REF] Fang | The Capacitated Lane Reservation Problem in Transportation Network[END_REF] extended a capacitated lane reservation problem. However, their works don't respect to the risk of transportation. In our work, the lane reservation strategy is applied to the hazardous materials transportation, where risk is the primary ingredient that separates hazardous materials transportation problems from other transportation problems.

It is well known that the likelihood of a traffic accident varies with road structure (e.g. the number of lanes and the use of lanes) [START_REF] Verter | A GIS-Based framework for hazardous materials transport risk assessment. Risk analysis[END_REF]. The risk of hazardous materials transportation though reservation lanes has much little impact on the population and environment. Nevertheless, lanes reservation will certainly impact the normal traffic because only special shipment can pass through the reserved lanes. Thus, decisionmaker (like a governmental planner, a regulator and/or manager of a shipper) should try to optimally catch a compromise between the level of risk and the impact on the normal traffic due to applying lane reservation strategy during the hazardous materials transportation. For the depictions above, it is understandable that the hazardous materials transportation problem based on lane reservation is a multiobjective optimization problem.

The contribution of this paper is that we address a quite new hazardous materials transportation problem based on lane reservation, which is to decide lanes to be reserved on the existing transportation network and select each hazardous materials shipment path from the reserved lanes. We first present a multi-objective integer linear programming model for the problem with multiple types of hazardous materials and multiple origin-destination pairs. We then apply the εconstraint method to solving the proposed problem. The remainder of the paper is organized as follows. In Section II, a formal problem description is given and the mathematical model is developed, respectively. In Section III the εconstraint method is adopted to solve the multi-objective model. Computational results on sets of randomly generated instances are reported to demonstrate the effectiveness of the proposed approach in Section IV. Finally, Section V gives the conclusion.

II. THE MODEL

A. Problem Description

The regional network for hazardous materials transportation, G = (V, A), consists of a number of nodes V = {1, …, N}. The nodes may be the origins generating a great quantity of hazardous materials, towns with lots of population or the destinations denoting disposal plants. Nodes are connected to each other by means of transportation routes A = {(i, j)}, i, j  V. An arc (i, j) denotes a transportation route from node i to node j. Kinds of hazardous materials must be carried from the origins O = {o 1 , o 2 , …, o No } to its corresponding destination nodes D = {d 1 , d 2 , …, d Nd }. We call the shipment of hazardous material beginning at origin node o w as shipment w.

The problem is to decide lanes to be reserved on the transportation network and select each hazardous materials shipment path from the reserved lanes to ensure that each shipment can be finished from its source to destination within a given travel time and the given risk threshold value of each reserved lane, meanwhile optimally come to a compromise between the total impact on the normal traffic and the total hazardous materials transportation risk.

B. Assumptions and Notation

The following assumptions are given:  Note that any transport network can be defined so as to satisfy that the probability of the hazardous material accident on an arc is constant and directly proportional to distance of the arc.  Note that each of the hazardous materials accidents happens independently.  It is assumed that there are at least two lanes on a road.

The following sets are given:

G = (V, A) directed transportation network O = {o 1 , o 2 , …, o No } origin nodes D = {d 1 , d 2 , …, d Nd } destination nodes
Let us define: TS ij : the travel time on a reserved lane of arc (i, j) K ij : the number of lanes on arc (i, j) C ij : the cost of a reserved lane on arc (i, j) T w : the deadline to accomplish shipment w PT ij : the threshold value of the probability of the hazardous material accidents on arc (i, j) w PS ij : the probability of hazardous material accident during the shipment w on a reserved lane of arc (i, j) E ij : the population exposure in the bandwidth along arc (i, j)

There are two decision variables below.

1 if there is a reserved lane on arc ( , ) and the shipment passes the arc

0 otherwise ij A w xw ij        1 if there is a reserved lane on arc ( , ) 0 otherwise ij A y ij      C. Formulation
The mathematical formulations for the various objectives and constraints can be constructed as follows. It is a multiobjective integer linear programming model.

Minimize 1 (, ) f Cy ij ij ij A    (1)
Minimize 2 (, )

ww f EP S x ij ij ij W ij A    (2)
Subject to 1, , , ( , )

w x wW o O o j A ww oj w j       (3) 1, , , ( , ) w
w x wW d D i d A w id w i       (4) , , , , ( ,
) ww x xw W i o d i j A ww ij ji jj        (5)
, ( , ) ,

w x yi j A w W ij ij    (6)
, ( , )

, ww PS x PT i j A w W ij ij ij w      (7)
, ( , ) , (, )

w TS x T i j A w W w ij ij ij A      (8)
{0,1}, ( , )

w x ij A ij   (9)
{0,1}, ( , )

y ij A ij   ( 1 0 )
We define two main objectives for the problem. The objective (1) is to minimize the total impact on the normal traffic. The impact of reserved lanes from node i to node j on the normal traffic can be translated into the lane reservation cost C ij , which is defined by

1 T ij K ij 
according to [START_REF] Wu | Heuristic for lane reservation problem in time constrained transportation[END_REF]. The objective (2) is to minimize the total hazardous materials transportation risk associated with a unit flow of each shipment. As point out by Erkut and Verter [START_REF] Erkut | Modeling of Transport risk for Hazardous Materials[END_REF], researchers have not come to a consensus on how to model the risk associated with hazardous materials transportation. In quantitative risk assessment, risk is defined as the product of the probability of hazardous materials accident and its consequence. One of the simple methods for risk is to that multiplying the hazardous materials accident probability by population exposure. Risk assessment is a quite important task but it is out of the aim of our work. Therefore, in this study, our proposed model uses this method as a surrogate for risk measure. How to formulate the risk for entire path between an origin and a destination can be found in [START_REF] Erkut | Modeling of Transport risk for Hazardous Materials[END_REF].

We introduce the following constraints which have to be satisfied to design an appropriate transportation network. Constraints (3) and ( 4) state that there is only one path for each shipment w starting from the origin node o w or arriving at the destination node d w . Constraint (5) ensures the flow conservation. Constraint [START_REF] Pesatori | Cancer incidence in the population exposed to dioxin after the "Seveso accident": twenty years of follow-up[END_REF] represents that there is a reserved lane in the path of one shipment if and only if this lane is reserved. Constraint [START_REF] Erkut | Chapter 9, Hazardous Materials Transportation[END_REF] specifies that the risk caused by all the shipment pass arc (i, j) can not exceed the risk threshold value of the corresponding arc. Constraint (8) guarantees that the total travel time for shipment w does not exceed the deadline T w . Constraints ( 9) and ( 10) specify that sign restrictions on the variables.

III. THE ΕPSILION -CONSTRAINT METHOD FOR THE PROPOSED MULTI-OBJECTIVE MODEL

Without loss of generality, multi-objective optimization can be formulated as follows:

min [f 1 (x), f 2 (x),..

., f k (x)] S

x 

where k is the number of objective function and S is the feasible solution space. A vector *

S

x  is said to be weakly Pareto optimal or a weakly efficient solution for the multiobjective problem if and only if there is no

S x  such that * () ( ) ii f xf x  for all   1,...,k i  . A vector * S
x  is said to be Pareto optimal or a strict efficient solution for the multiobjective problem if and only if there is no

S x  such that * () ( ) ii f xf x  for all   1,...,k i 
and at least one strict inequality. The image of all the efficient solutions is called Pareto front. One of the multi-objective optimization techniques is the εconstraint method (also called the e-constraint or trade-off method) introduced by Haimes et al. [START_REF] Haimes | On a bicriterion formulation of the problems of integrated system identification and system optimization[END_REF]. This method is based on minimizing of one of the objective functions and restricting the other objectives within some allowable values ε i . Hence, a single objective minimization is carried out for the most relevant objective function subject to additional constraints on the remaining objective functions. It generates one point of the Pareto front at a time, and then the levels ε i are altered to generate the entire Pareto front. Miettinen in 1999 [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF] proved that if it exists, a solution to the ε-constraint formulation is weakly Pareto optimal.

As Section II presented, the proposed multi-objective model for hazardous materials transportation problem based on lane reservation has two competing objective functions subject to a set of constraints. The first objective function is the cost of the lane reservation within the transportation network, while the second objective function is the total transportation risk on condition that those lanes that all the shipments passed are reserved. That is to say, a reserved lane is a condition of a shipment. Therefore, we can consider the first objective function as the main objective function. Hence, the proposed multi-objective problem can be transformed into the single objective problem with the objective function 1 f being minimized while the objective function 2 f being imposed as a constraint to the feasible solution space. Fig. 1 depicts the flowchart of the proposed ε-constraint method, which is explained below.

1) The problem presented in Section II (called Problem P) is transformed into Problem P 0 . In Problem P 0 , the objective is to minimize the 1-th objective function of P, and the 2-th objective function of P is dealt with as a constraint of P 0 . P 0 : Minimize (, )

Cy ij ij ij A   (1)
Subject to , (, )

ww EP S x ij ij ij W ij A     (11)
and constraints ( 3)-( 10). Where ε is an upper bound of the value of the objective function 2 f and is not necessarily a small value close to zero. 2) A systematic variation of ε yields a set of weakly Pareto optimal solutions. However, improper selection of ε can result in a formulation with no feasible solution. Note that ε is bounded by [lower limit, upper limit]. We define the limits according to a) and b). a) for the objective function 1 f , an optimal vector

( * ij y , * w ij
x ) can be found by solving the following Problem P 1 .

P 1 : Minimize (, ) Cy ij ij ij A   (1)
Subject to constraints (3)- [START_REF] Erkut | A bi-level programming application to hazardous material transportation network design[END_REF].

ww EP S x ij ij ij W ij A   (2)
Subject to constraints (3)- [START_REF] Erkut | A bi-level programming application to hazardous material transportation network design[END_REF].

And let 2 ' ()

w ij f x be the lower limit of ε, i.e., 2 '
()

w ij f x   .
3) Choose randomly and uniformly  between

IV. NUMERICAL RESULTS

The proposed model formulation in Section III was implemented in C. All the computations were done on a HP PC with a Pentium IV Processor 3.0 GHZ. The mixed integer programming (MIP) solver provided by the commercial optimization software package CPLEX (Version 12.1) was applied to solving the integer programming (IP) problems, i.e., Problems P 0 , P 1 and P 2 . CPLEX was run in default setting.

A. Experimental Setup

In order to test the proposed model, we design the random test problems as follows. The transportation network was generated based on the random network topology generator introduced by Waxman [START_REF] Waxman | Routing of multipoint connections[END_REF]. In his model, given the number of the nodes, the nodes of the network are randomly and uniformly distributed in the plane and the arcs are added according to probabilities that depend on the distances between the nodes. The probability function to have an arc between two nodes can be found in [START_REF] Waxman | Routing of multipoint connections[END_REF]. We generate the origin-destination (o-d) pairs randomly from the nodes of the network.

All the test instances were generated with a random number stream. Let Uniform (a, b) stand for a uniform distribution between parameters a and b. We generated several sets of randomly generated test instances using the parameters described below. T i,j and w PS ij respectively denotes the travel times and the probability of accident of the hazardous material shipment w on arc (i, j) without reserved lane. T ij were generated by using the following parameters:

T i,i+1 = T i+1,i = Uniform (2, 4), T i,j =T j,i = 1 ,1 j T kk ki    
and the travel times TS ij = T i,j * Uniform (0.5, 0.8). The deadline T w is set to be distance and road condition on the accident probability given by [START_REF] Qiao | Utilization of accident databases and fuzzy sets to estimate frequency of HazMat transport accidents[END_REF].

w P ij = w P ji = 1 ,1 j w P kk ki     and w PS ij = w P ij * Uniform (0.2, 0.5),
where one unit of risk is 10 -7 events * yr -1 . Note that PT ij = W w P ij w  * Uniform (0.4, 0.8). Lastly, the population exposure in the bandwidth along arc (i, j), E ij , is generated by Uniform (10, 80), whose unit is 10 4 ; K ij , the number of lanes on arc (i, j), can be defined by Uniform [START_REF] Leonelli | New detailed numerical procedures for calculating risk measures in hazardous materials transportation[END_REF][START_REF] Erkut | A bi-level programming application to hazardous material transportation network design[END_REF] in accordance with [START_REF] Qiao | Utilization of accident databases and fuzzy sets to estimate frequency of HazMat transport accidents[END_REF].

B. Experiment for Randomly Generated Instances

Firstly, let us define the network density as D =

(1 )

A VV  , 0  D  1,
where |V| and |A| is the number of nodes and arcs respectively [START_REF] Wu | Heuristic for lane reservation problem in time constrained transportation[END_REF]. The maximum number of arcs is |V|(|V|-1), so the maximal and minimal density is 1 (for complete networks) and 0 respectively. For each problem, we allowed the software to run until it was solved to optimality.

Let n = 20. Hence, for each test, 20 different values of  were generated so that 20 Pareto optimal solutions could be obtained each time. Total computational time means the computational time for obtaining 20 solutions in all, and average computational time is defined as total computational time divided by 20. The detailed computational results are analyzed as follows.

Tables I and II respectively gives the summary result on the random generated instances with D  0.35 and 0.35< D  0.55. As can be seen from Tables I andII, increases in the number of the nodes, arcs and shipments of the problems lanes led to a steady increase in the number of binary variables, the number of constraints and the average computational time (in CPU seconds) for finding a Pareto optimal solution for the instances. In more detail, both in Tables I and II, given the number of the shipments, the larger the number of nodes the more time spent on finding an optimal solution. Table Ⅰ represents that given the number of shipment, the time for set 1-5, 6-10, 11-15 and 16-20 raises gradually. Given the number of the nodes, the larger the number of shipment, the more time. In tableⅠthe average time for instance sets 7-10 is more than the counterpart for sets 11-14, respectively. Given the number of nodes and shipments, the larger the number of arcs, the more time. Take instances sets 6-10 with nodes from 30 to 70 and 10 shipments as an example, the average computational time for the bigger network density is more than the counterpart for the smaller one. The number of binary variables, the number of constraints work as the time in the same way, and analyses on them are omitted here. Then, we also have confirmed whether the solutions of a problem could be effective by taking instance set 7 and set 37 as an example. For this test, 20 different values of  were generated, and 13 and 19 different Pareto optimal solutions could be obtained in all, as shown in Fig. 2 and Fig. 3. As these figures indicate, all the solutions seem likely to make up of a curve. Decision-makers can choose one of points as a Pareto optimal solution for the problem according to their preference.

V. CONCLUSIONS

In this paper, we presented a new problem based on the lane reservation strategy in the hazardous materials transportation network. We developed a multi-objective integer programming model for the problem and then apply the ε-constraint method to solving the proposed model. Numerical results of a large set of examples showed the effectiveness of the ε-constraint method. The ε-constraint method can solve the proposed problems within reasonable time and moreover generate candidate solutions that can shape an efficient Pareto front for the multi-objective problem.

Although our problem could be solved optimally by CPLEX, we may have to resort to other algorithms, with which it may be possible to solve even larger problems within shorter time. This may be a future research direction. Another future work may be to extend the proposed mathematical formulation. For example, various objectives of the hazardous materials transportation problem can be implemented in our model with consideration of the risk equity.

  be the upper limit of ε, i.e., be found by solving the following Problem P 2 .

4 )

 4 . Solve Problem P 0 , obtain a Pareto optimal solutions. Repeat 3) n times by varying  Pareto optimal solutions, which can theoretically form a Pareto front.

Figure 1 .

 1 Figure 1. The flow chart of the ε-constraint method

(

  o w , d w ) * Uniform ( 2 , 3 ), where distance (o w , d w ) is the shortest travel time from o w to d w in a reserved path. Similarly, let ,1 w P ii = 1, w P ii  = Uniform (8, 20), considering the effects of the number of lanes, truck configuration, population density, Solving Problem P 1 and obtaining the upper limit on ε Solving Problem P 2 and obtaining the lower limit on ε Choosing a value of ε Solving problem P 0 and obtaining a Pareto optimal solution Obtaining a Pareto Front n times

Figure 2 .Figure 3 .

 23 Figure 2. Solutions of the problem with 40 nodes, 10 shipments and D 0.35

TABLE I THE

 I

	T ABLE II THE SUMMARY RESULT ON THE RANDOM GENERATED INSTANCES
	WITH 0.35< D 0.55
	SUMMARY RESULT ON THE RANDOM GENERATED INSTANCES
	WITH D 0.35

This work was partially supported by the Cai Yuanpei Program between the French Ministries of Foreign and European Affairs and the Higher Education and Research and the Chinese Ministry of Education, under Grant No. 24021SH.