F Cazals

T Dreyfus
email: tom.dreyfus@inria.fr

S Sachdeva
email: sachdeva@cs.princeton.edu

N Shah

Greedy Geometric Optimization Algorithms for Collection of Balls

Keywords: Geometric optimization, geometric approximation, collection of balls, inner approximation, molecular conformations Optimisation géométrique, approximation géométrique, collections de boules, approximation par défaut, conformations moléculaires

Modeling 3D objects with balls is routine for two reasons: on the one hand, the medial axis transform allows representing a solid object as a union of medial balls; on the other hand, selected shapes, and molecules in particular, are naturally represented by collections of balls. Yet, the problem of choosing which balls are best suited to approximate a given shape is a non trivial one. This paper addresses two problems in this realm. The rst one, conformational diversity selection, consists of choosing k molecular conformations amidst n, so as to maximize the geometric diversity of the k conformers. The second one, inner approximation, consists of approximating a molecule of n balls with k < n balls. On the theoretical side, we demonstrate that for both problems, a geometric generalization of max k-cover applies, with weights depending on the cells of a surface or volumetric arrangement. Tackling these problems with greedy strategies, it is shown that the 1 -1/e bound known in combinatorial optimization applies in some cases but not all. On the applied side, we present a robust and eective implementation of the greedy algorithm for the inner approximation problem, which incorporates the calculation of the exact Delaunay triangulation of a points whose coordinates are degree two algebraic number, of the medial axis of a union of balls, and of a certied estimate of the volume of a union of balls. In particular, we show that the inner approximation of complex molecules yields accurate coarse-grain models with a number of balls 100 times smaller than the number of atoms, a key requirement to simulate crowded protein environments.

Modeling complex 3D shapes is commonplace in science and engineering, and simple primitives such as balls play a central role in this process, for two reasons. On the one hand, the medial axis transform (MAT) allows representing a shape as a collection of balls [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], usually innite, so that sub-sampling such balls naturally yield approximations. On the other hand, (hierarchical) models represented by balls are ubiquitous, for example in molecular modeling, but also in robotics, computer graphics and CAGD, where bounding sphere hierarchies provide an elegant way to perform fast and numerically reliable collision detection. In this context, this paper addresses the following two problems, which we phrase using the molecular modeling terminology, even though the semantics of the balls might be dierent:

Conformational selection. Given a set of n conformations of a molecule represented by a collection of balls, select a diverse ensemble consisting of s < n conformations.

Inner approximation. Given a (molecular) model consisting of n balls, provide an accurate volumetric approximation of this model using s < n balls, contained in the original model.

These two problems are actually connected to a variety of research veins, namely (i) geometric approximation algorithms for 3D shapes, (ii) medial axis constructions and Voronoi diagrams, (iii) (geometric) approximation algorithms in general and max k-cover in particular, (iv) robust geometric software development, and (v) applications in structural biology. We now briey comment on recent work in these directions.

As already mentioned, the problem of approximating 3D shapes is related to the medial axis transform (MAT). The particular case of a shape bounded by a smooth surface motivated the introduction of the MAT approximation using medial balls centered on specic Voronoi vertices called poles [START_REF] Amenta | Accurate and ecient unions of balls[END_REF], an idea later re-used to approximate a shape bounded by a triangulated surface [AAK + 09, SKS12]. This MAT approximation was also used for the sphere-tree construction [START_REF] Bradshaw | Adaptive medial-axis approximation for sphere-tree construction[END_REF], a representation to perform hierarchical object modeling and collision detection, and to improve the grasping quality in robotics [START_REF] Przybylski | Unions of balls for shape approximation in robot grasping[END_REF]. For a shape with smooth boundary, the previous MAT approximation typically comes with a guarantee, namely the convergence of the Hausdor distance between the input boundary and that of the approximation. In a broader context, the problem of approximating a bounded open set has also been investigated recently. In [START_REF] Giesen | The scale axis transform[END_REF] the authors introduce the scale-axis transform, which consists of scaling forth and back medial balls, so as to simplify a shape representation. It is worth noticing that all the works just mentioned rely on Voronoi diagrams, generally for the Euclidean distance, but also for a multiplicative distance in [START_REF] Giesen | The scale axis transform[END_REF]. Consequently and from an implementation perspective, geometric algorithms from the Computational Geometry Algorithms Library [cga] (CGAL), but also number types from the LEDA [START_REF] Mehlhorn | LEDA: a platform for combinatorial and geometric computing[END_REF] and CORE [START_REF] Karamcheti | A core library for robust numeric and geometric computation[END_REF] libraries play a key role to implement such algorithms.

Our problems are also related to approximation algorithms in general, and greedy strategies in particular. As we shall see, of particular interest for both problems is max k-cover, which cannot be approximated within a ratio of 1 -1/e + ε unless P = NP [START_REF] Feige | A threshold of ln n for approximating set cover[END_REF].

Last but not least, our incentive to tackle these two problems comes from computational structural biology, whose ultimate goal is to unravel the relationship between the structure and the function of macro-molecules and macro-molecular machines. Originating with the work of Richards [START_REF] Lee | The interpretation of protein structure: Estimation of static accessibility[END_REF], molecular models represented as collections of van der Waals (vdW) balls and associated ane Voronoi diagrams have been instrumental to describe atomic packing properties [MJLC87, MLJ + 87], to compute and decorate molecular surfaces [START_REF] Connolly | Analytical molecular surface calculation[END_REF][START_REF] Akkiraju | Triangulating the surface of a molecule[END_REF], to ex-Inria hibit correlations between structural and biological -biophysical properties of protein interfaces [START_REF] Bahadur | Dissecting subunit interfaces in homodimeric proteins[END_REF][START_REF] Malod-Dognin | Characterizing the morphology of protein binding patches[END_REF], to select diverse conformational ensembles for mean eld theory based docking algorithms [LSB + 11], or to nd entrance / exit passages to active sites [YFW + 08]. More recently, compoundly weighted Voronoi diagram [START_REF] Cazals | Multi-scale geometric modeling of ambiguous shapes with toleranced balls and compoundly weighted α-shapes[END_REF][START_REF] Dreyfus | Assessing the reconstruction of macro-molecular assemblies with toleranced models[END_REF] proved instrumental to assess the reconstruction of molecular machines involving up to circa 450 polypeptide chains [ADV + 07], these reconstructions being plagued with uncertainties on the shapes and positions of the proteins. In fact, the inner covering problem is directly related to the design of simplied yet accurate protein models, in particular in the perspective of simulating crowded whole cellular environments [START_REF] Mcguee | Diusion, crowding and protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF][START_REF] Goodsell | The machinery of life[END_REF] see also the beautiful illustrations of D. Goodsell1 .

Contributions and Paper Overview

Our two problems are concerned with shapes represented as a union of balls. These balls may have a molecular semantics, or may be carry a purely geometric meaning, e.g. in the MAT of any 3D shape. In using balls, a key advantage is that the structure of the medial axis of the union of these balls is known exactly [AM97, AK01]: it is actually coded by the α-shape of the input balls [START_REF] Edelsbrunner | Weighted alpha shapes[END_REF], and the Voronoi diagram of the points found on the boundary of the union. In this context, our contributions are threefold. On the theoretical side, we demonstrate that for both problems, a geometric generalization of max k-cover applies, with weights depending on the cells of a surface or volumetric arrangement. In particular, it is shown that the 1 -1/e bound known in combinatorial optimization applies in some casesbut not all. From a geometric approximation perspective, our results depart from previous work since we focus on an approximation guarantees obtained with a nite set of balls rather than asymptotically. Second, on the software development side, we present an eective implementation for the inner approximation problem, based on state-of-the art geometric software. Finally, for the inner approximation problem, we present accurate coarse-grain protein models using a number of balls 100 times smaller than the number of atoms.

Geometric Optimization Problems for Balls

In the sequel, a molecule stands for a collection of balls. For the conformational selection, the molecule is termed a conformer when the relative position of the balls may change. The generic term conformation refers to a molecule or to a conformer, and we shall manipulate a collection of conformations C = {C i } i=1,...,n . The 3D domain spanned by the conformations is the union of their dening balls, and is denoted F C = ∪ Ci∈C C i . For the inner approximation, we consider a single molecule. The set C refers to its constituting balls, and the domain F C is the union of these balls.

Geometric pre-requisites

Volume and surface decompositions.

The spheres bounding the balls of a collection of conformations induce two decompositions: a decomposition of the volume occupied by the balls; and a decomposition of each sphere into spherical patches. More precisely, the decomposition of the volume F C induced by the spheres is called a volumetric arrangement (or volumetric decomposition). This arrangement consists of a collection of cells A = {A i } such that the interior of each cell is connected. Each such cell is bounded by 2D cells, called surface patches, found on the spheres bounding the balls. On a given sphere, these patches are induced by the intersection circles with neighboring spheres. The collection P = {P i } of all such patches denes a surface arrangement (or surface decomposition). A surface patch of a sphere which is not contained in any other ball is called exposed, and so is its supporting sphere. The collection of all such patches makes up the boundary ∂F C of the domain F C . Note that each patch is bounded by circle arcs which are themselves delimited by points (generically) found at the intersection of three spheres. See Fig. 1 for a 2D illustration.

Medial axis of a collection of balls.

The medial axis of the domain F C consists of the points having several neighbors on ∂F C . Let a boundary point of F C be a 0-cell of ∂F C . As proved in [START_REF] Amenta | The medial axis of a union of balls[END_REF] and illustrated on Fig. 2, the MA consists of so-called singular simplices of the α-complex for α = 0, together with a subset of the Voronoi diagram of the boundary points located within the α-shape.

Generic Geometric Optimization Problems

Conformational selection. We shall be concerned with two classes of combinatorial optimization problems. To state them from a combinatorial viewpoint (see section 3.1 and 3.2 for the connexion with conformations), assume we are given a base set U = {U i } i=1,...,m of interior disjoint cells (think cells of the volume or surface arrangement), and a collection of sets C = {C i } i=1,...,n called the pool (think conformations), where each set is a union of cells. For a subset S ⊂ C, denote F S = ∪ Cj ∈S C j the union of the sets in S. Cells and sets shall be subsets of R 3 , so that the inclusion of a cell U i in a set C j is naturally dened.

For the rst class of problem, assume we are given a weight function w, i.e. a real valued function dened over the cells. Let C s stand for the subsets of C of size s. We dene:

Problem 1. Given a weight function w, nd a subset Ŝ of C of size s, called the selection, such that:

Ŝ = arg max S∈( C s )
w(S), with w(S) = Ui⊂∪ S Cj w(U i ).

(1)

For the second class of problems, assume the weight function depends not only on the cells of the decomposition, but also on the selection S, which we denote w S (U i ). We wish to solve: Problem 2. Given a weight function w S , nd a subset Ŝ of C of size s, called the selection, such that:

Ŝ = arg max S∈( C s ) w(S), with w(S) = Ui⊂∪ S Cj w S (U i ).
(2) Inner approximation.

Consider a collection of n balls making up the domain F C , that is F C = ∪ Bi∈C B i , and let Vol(D) denote the 3D volume of a domain D. We introduce the following covering problem, phrased as an inner approximation problem: Problem 3. Find a collection S of s balls, such that F S ⊂ F C and Vol(F C \F S ) is minimized. This problem is more constrained than the minimization of the volume of the symmetric dierence Vol(F S ∆F C ). Yet, the inner approximation problem is natural since it requires using balls centered on the medial axis of the input domain, since any other ball would be contained within a maximal ball. Inria Complexity issues.

Our problems are intimately related to max k-cover. Given a set U of n points, and a collection C of subsets of U, max k-cover aims at selecting k subsets from C so as to maximize the number of points from U which get covered [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF][START_REF] Feige | A threshold of ln n for approximating set cover[END_REF]. (In the literature, this problem is sometimes called set cover [START_REF] Pc Fishburn | Pick-and choose heuristics for partial set covering[END_REF]. To avoid confusion, we consider that the set cover problem aims at minimizing the number of sets in C to cover at least k elements from U.) We note that the classical max k-cover is a special case of Problem 1 with function w assigning a unit weight to all cells. Since max k-cover is a NP complete problem, a polynomial time solution both |C| and s cannot be expected. However, the problem is in P for a xed s since all subsets of size s can be probed. But this brute force method is doomed to fail even for moderate s, which calls for alternate strategies, the greedy strategy being the most natural one.

2.3

The Greedy Strategy

The greedy strategy consists of s iterations, the kth step consisting of selecting the C j maximizing the weight of the union of the C j . The selection achieved at step k shall be denoted S k .

Unfortunately, the selection S s may not realize the optimum solution, as illustrated by the simple example of Fig. 3. Thus, the performance assessment of greedy relies on the worst-case ratio between the solution returned and the optimal one. For max k-cover, this ratio is known to be of 1 -1/e, and is tight [CFN77, NWF78, FG89, Fei98].

Conformational Selection and the Greedy Strategy

We rst examine problems concerned with the selection of a subset of conformers maximizing some diversity criterion.

Volumetric Decompositions and the Greedy Strategy

Problem 1 from volumetric decomposition.

Consider the base set A whose cells are those of the 3D arrangement. In Eq. (1), let w be some general function dened on the cells of the volumetric decomposition. The weighting scheme is called non-negative provided all weights are ≥ 0the most natural example being the standard Euclidean volume.

The approximation ratio of the greedy strategy and its optimality are usually proved in the uniform weight case [START_REF] Cornuejols | Location of bank accounts to optimize oat: An analytic study of exact and approximate algorithms[END_REF][START_REF] Nemhauser | An analysis of approximations for maximizing submodular set functions[END_REF][START_REF] Pc Fishburn | Pick-and choose heuristics for partial set covering[END_REF][START_REF] Feige | A threshold of ln n for approximating set cover[END_REF]. We shall prove the following Theorem 3.1. Consider a volumetric decomposition with non-negative weights. For Problem 1, the greedy approach has an approximation ratio of 1

-(1 -1/s) s > 1 -1/e.
We shall use the following notation. The conformer selected at the k th step is denoted C k , and the weight of the optimum set of conformers OP T . Also, let w * (C k ) be the sum of the weights of the new elements in C k that have not been covered in C j , 1 ≤ j < k (i.e. the weight increment at step k). We start with a lemma needed to prove theorem 3.1. The proof of the lemma follows the usual one for max k-cover, but we include it for two reasons: rst, it helps spotting the condition on the weight-functions w (the positivity is mandatory); second, in section 4.2, we shall re-use the skeleton of this lemma to prove a result of the greedy strategy for inner covering, with respect to the total volume of F C instead of the optimum. Lemma 3.2. For 1 ≤ k ≤ s, the following holds:

w * (C k ) + 1 s k-1 j=1 w * (C j ) ≥ OP T s .
(3)

RR n°8205

Proof. At the k th step, we select C k that maximizes the weight of the new cells U i being covered. Because the cells selected up to step k -1 may cover cells which are not covered by OPT, the weight of the cells that are covered by the optimum solution but not yet covered by the (k -1)

is at least OP T - k-1 j=1 w * (C j ) (4)
Since w is non-negative, the union-bound property states that for any collection of conformers C 1 , . . . , C p , one has w(C

1 ∪ • • • ∪ C p ) ≤ i=1,...,p w(C i ).
Since all the cells involved in Eq. ( 4) are covered by the optimum set of conformers, by the union-bound property, there must exist one conformer,not yet selected, that covers these new cells with total weight at least

1 s   OP T - k-1 j=1 w * (C j )   .
(5)

Since C k maximizes the weight of the new cells being covered, we must have

w * (C k ) ≥ 1 s   OP T - k-1 j=1 w * (C j )   . (6)
Rearranging completes the claim.

Remark 1. The non-negativity assumption is critical in the proof of Lemma 3.2. As a counterexample, consider the sets C 1 = {e1, e2}, C 2 = {e2, e3} with w(e1) = w(e3) = 1 and w(e2) = -1.

The union-bound fails for w(C 1 ∪ C 2 ). This remark is of particular interest in bio-physics, where atoms are decorated with physical, chemical or biological properties. For example, a weighting function that would take into account the electrostatic properties, which may be negative, would preclude the application of the previous lemma.

Using Lemma 3.2, the proof of Thm. 3.1 goes as follows:

Proof. Multiplying the inequality obtained in the previous lemma by (s -1)/s and adding to the inequality for step two, we get

w * (C 1 ) + w * (C 2 ) ≥ 1 + s -1 s

OP T s

We multiply this equation again by s-1 s and add to the equation for step three, and so on. We get the following,

k j=1 w * (C j ) ≥ 1 - s -1 s k OP T For k = s, we get, s j=1 w * (C j ) OP T ≥ 1 - s -1 s s
The left hand side is the ratio of the weight of the subset of C chosen by the greedy approach and the optimum solution i.e. that approximation factor and hence we have the above theorem.

The fact that the above ratio is greater than 1 -1 e for all s is a trivial exercise.

Inria

We now prove that the bound of the previous is tight (see also Fig. 4):

Theorem 3.3. The greedy approach cannot perform better than 1 -(1 -1/s) s .

Proof. Fix a given s. We shall construct an example where the greedy approach can achieve an approximation ratio arbitrarily close to

1 -(1 -1 s ) s . Let A = {A i } i=1,...,(s 2 +s) ∀i, j s.t. 0 ≤ i < s, 1 ≤ j ≤ s, w(A i.s+j ) = 1 s 2 s -1 s i ∀j s.t. 1 < j ≤ s, w(A s 2 +j ) = 1 s s -1 s s -ǫ
The conformers are dened as follows

C = {C i } i=1,...,2s ∀i s.t. 1 ≤ i ≤ s, C i = i.s j=(i-1).s+1 A j ∀i s.t. ≤ i ≤ s, C s+i = j≡i (mods) A j
Simple calculations lead us the following total weights for the conformers

∀1 ≤ i ≤ s, w(C i ) = 1 s s -1 s i-1 ∀1 ≤ i ≤ s, w(C s+i ) = 1 s -ǫ
The optimum choice of S with |S| = s is clearly {C i } i=s+1,...,2s with total weight 1 -sǫ, whereas the greedy method would choose {C i } i=1,...,s , with a maximum weight of 1 -(1 -1 s ) s , giving an approximation factor is arbitrarily close to 1 -(1 -1 s ) s .

Surface Decompositions and the Greedy Strategy

Problem 2 from surface decomposition. Consider the base set P = {P i } whose cells are those of the 2D arrangements induced on each sphere by the intersection circles with all the other spheres, computed e.g. with the algorithm from [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF]. Special cells of this arrangement are those which are exposed, i.e. contribute to the boundary of the union of balls. Focusing on these patches yields an instantiation of Problem 2, the dependence upon the selection S consisting of discarding the patches which are not exposed with respect to the selection. For example, w S (P i ) = surface area of patch P i i P i is found on the boundary of the union F S , and 0 otherwise. Interestingly, maximizing the boundary surface of the selection is an indirect way to ascertain some diversity, since the overlap between conformers is minimized. Notice, though, that as opposed to the volume, the boundary surface area is not a monotonic function of the number of conformers. That is, for two selections S 1 and S 2 with S 1 ⊂ S 2 , one has volume(S 2 ) ≥ volume(S 1 ), a property that may not hold for the boundary surface area. Surface decomposition, boundary surface weight w S . For volumetric decompositions, the previous bound indicates that one is always above 63% (1 -1/e) from the optimum. Unfortunately, such a result does not hold for problem 2: RR n°8205

Observation 1. Consider a surface decomposition. For Problem 2, the greedy approach may have a worst-case approximation ratio as bad as 1/s 2 .

Proof. Consider a large ball B, and place s small non-intersecting balls (B 1 , . . . , B s ) with their centers on the surface of B. The surface of each B i is now divided into 2 patches. To the patch which lies inside B, we assign a weight of s. To each surface patch of B covered by some B i , we assign a weight of 1 + ǫ. All other surface patches are assigned a weight of 0.

The greedy strategy would rst pick B because it has the largest exposed weight of s(1 + ǫ). Now picking any s -1 of the B i 's would leave us with an exposed weight of only s(1 + ǫ) -(s -1)(1 + ǫ) = 1 + ǫ. On the opposite, a selection of the s smalls balls would have given us total exposed surface weight of s 2 . This approximation factor is arbitrarily close to 1/s 2 .

Inner Approximation and the Greedy Strategy

We establish two results for the inner approximation: rst the performance of greedy with respect to the optimal solution, and second, with respect to the total volume Vol(F C ).

Approximation Factor

As discussed when introducing problem 3, the inner approximation problem requires using balls centered on the medial axis. But the medial axis is a cell complex with two dimensional faces, so that one has an innite collection of balls to choose from. To circumvent this diculty, consider the following classical lemma [Ber87]:

Lemma 4.1. Consider two intersecting spheres Σ 1 and Σ 2 in 3D, and dene their convex linear combination, namely

Σ λ = λΣ 1 + (1 -λ)Σ 2 , with λ ∈ [0, 1].
The ball bounded by Σ λ is contained in the union of the balls bounded by Σ 1 and Σ 2 .

Denote B * p a maximal ball centered on a vertex p of the medial axis, and let V be the set of vertices of the medial axis of F C . By the structure theorem of the medial axis of a union of balls [START_REF] Amenta | The medial axis of a union of balls[END_REF], this set is nite. We shall use this set to run the greedy algorithm, based on the following:

Observation 2. The input domain F C satises

F C = i B i = v∈V B * v . (7)
Proof. We shall prove that any maximal ball B * p is contained in the union of at most three balls centered on vertices from V. Omitting the trivial case of a singular vertex of the medial axis, we rst note that there are three cases to be analyzed, namely when p belongs to a singular edge of the medial axis, when it belongs to a singular triangle, or when it belongs to a (possibly clipped) Voronoi face f . Case 1. This is exactly the case covered by lemma 4.1. In this case, the portion of the pencil contains the intersection circle between the two spheres dening the singular edge. Case 2. The second case contains two sub-cases, namely when p lies in the interior of a Voronoi edge, and when p lies in the interior of the Voronoi facet f . The rst sub-case is again the case of lemma 4.1 all the spheres in the portion of the pencil contain the three boundary points dening the Delaunay triangle dual of the Voronoi edge in question. For the second one: let c be Inria any Voronoi vertex of f belonging to V, let L be the ray emanating from c and passing through p, and let d be the intersection point between L and the boundary ∂f of f . Point d belongs to either a Voronoi edge or to an α-shape edge (if the Voronoi facet is a clipped Voronoi facet in the medial axis). Call e and f the endpoints of this 1-cell of the medial axis. Now, by lemma Case 3. Amenable to the analysis carried out for Case 2.

Thus, since any maximal ball is contained in the union of at most three balls centered at vertices from V, the claim holds.

The corollary of the previous observation is that the balls centered on the medial axis constitute the pool of candidate balls to choose from by the greedy strategy. Thus, theorem 3.1 applies, that is:

Theorem 4.2. For Problem 3, the greedy approach based on the maximal balls centered on the vertices of the medial axis has an approximation ratio of 1 -(1 -1/s) s > 1 -1/e.

Worst Case Bound with respect to the Total Volume

Approximation bound.

We generalize the previous result with respect to the total volume of the input domain:

Lemma 4.3. Let V = Vol(F C ) be the volume of union of given n balls and let GREEDY be the volume of union of s balls selected by the greedy algorithm. These volumes satisfy:

GREEDY V ≥ 1 -1 - 1 n s (8)
Proof. In the proof of the approximation factor of greedy algorithm for volumetric decomposition given in lemma 3.3, note that it is valid for any solution and not only the optimum solution, i.e.

no property of the optimum solution is required. Thus we replace the optimum solution by a solution selecting the given n balls. Thus we get the following equation.

w * (C k ) ≥ 1 n V - k-1 i=1 w * (C i )
where C i is the i th ball selected by the greedy algorithm, and w * (C k ) is the new volume of C k not covered by any of C i , 1 ≤ i < k. Solving it in the manner similar to that used in the proof of Them 3.3 yields:

GREEDY = s i=1 w * (C i ) ≥ V • 1 -1 - 1 n s Tight example.
Consider n balls of same radii. Then greedy algorithm would select any s balls out of it. This would contribute a volume equal to s/n times the total volume. Also note that

s n = 1 -1 - s n ≈ 1 -1 - 1 n s
for very large values of n. In fact, this is the best that can be done, i.e. no algorithm can approximate union of n balls with approximation factor greater than s/n in worst case.

RR n°8205 5 Inner Approximation: Software and Experiments

We focus on the inner approximation since the implementation of the conformational selection merely requires a robust algorithm to compute the volume of a union of balls in the volume case, and an algorithm to compute arrangements of all types of circles on a sphere in the surface case.

Greedy Algorithm

Basic algorithm.

The input consists of a collection of balls, and of a selection size s or a target ratio τ (the volume of the selection divided by the volume of the input balls should be at least τ ). The output consists of an ordering of the selected balls, together with the increment in volume associated to each ball. We also report the Betti numbers along the selection, computed with the algorithm from [START_REF] Delnado | An incremental algorithm for Betti numbers of simplicial complexes on the 3-sphere[END_REF].

The algorithm consists of iteratively selecting the ball providing the best volume increment, selected from a priority queue containing all candidates. Upon selecting ball say B i , we recompute the volume increments of all candidate balls intersecting B i .

Imposing connectedness of the selection.

For selected applications, the domain F S should be connected: for example, the selection associated to a connected molecule should also be connected. To meet this constraint, the following heuristic is used. Let S k be the selection upon termination, and consider the exposed balls i.e. the balls contributing to the boundary ∂F S . Split these balls into two groups L and L c , namely the largest component (in number of exposed balls), and the remaining ones. We aim at connecting L to one of the connected components of L c . To do so, using the Delaunay triangulation of the centers of the balls in L ∪ L c ∪ (S\S k ), we compute the shortest path joining a center of a ball in L to a center of a ball in L c . This shortest path uses centers of balls in S\S k , which are added to the section. This process is iterated until one connected component remains.

Geometric Objects

The previous algorithm involves elaborate geometric objects, which we present now by following the ow of the algorithm, mentioning the CGAL2 classes used and their template parameters when appropriate.

The Delaunay triangulation DT B of the input balls, and the associated α-shape. Following classical usage, we call K the kernel used to instantiate the CGAL classes Regular_triangulation_3 and Alpha_shape_3. Two options for K are discussed below.

The Delaunay triangulation DT V of the boundary points of ∂F S . Two diculties are faced to construct DT V . First, more than three co-planar points are generic in DT V [START_REF] Amenta | The medial axis of a union of balls[END_REF]. Second, since a boundary point is found at the intersection of three input spheres, its coordinates are degree two algebraic numbers. We therefore store these points using the CGAL spherical kernel Spherical_kernel_3 [START_REF] De Castro | Design of the cgal spherical kernel and application to arrangements of circles on a sphere[END_REF], instantiated with K. The two options for K, referred to as the inexact and the exact kernels in the sequel, are:

Exact_predicates_inexact_constructions_kernel, the underlying number type (NT) to store the coordinates of the boundary points being a double.

Exact_predicates_exact_constructions_kernel_with_sqrt, the underlying number type to store the coordinates being either CORE::Expr or LEDA::real. Additionally, a map is used to associate a singular or regular facet from the α-shape of DT B to each boundary point. To handle these diculties, we implemented a dedicated kernel denoted DTV_kernel, dening a new point type for the boundary points. This kernel is actually templated by two parameters:

First, a ball identier, used to record the three input spheres dening a boundary point. These identiers are used to handle the aforementioned special cases, so as to avoid the numerical calculation of a predicate whose sign can be inferred from the fact that the input points lie on a set of known input spheres. Practically and since an input ball corresponds to a vertex of the α-shape of DT B, the vertex handle of the α-shape is taken as identier.

Second, a number type used to represent the coordinates of the boundary points, the two options being the NT associated to the aforementioned inexact and exact kernels. One comment is in order about the Voronoi diagram DT V * dual of DT V , since medial ball associated to selected Voronoi vertices are used by greedy. With the inexact kernel, the input points of DT V are approximations of the exact boundary points, since the degree two algebraic number get converted to doubles. For these points, the combinatorial structures of DT V and DT V * are exact (exact predicates are used), but the embedding of the Voronoi vertices of DT V * is inexact (inexact constructions are used). With the exact kernel, the input points of DT V are exactly the boundary points. Moreover, the embedding of the Voronoi vertices is exact (exact constructions are used).

The medial-axis of the union of input balls. We store the medial axis as a container of polygons, possibly degenerate for singular vertices and edges of the α-shape [START_REF] Amenta | The medial axis of a union of balls[END_REF]. Our polygon class inherits inherit from the CGAL class Polygon_2 (embedded in 3D), instantiated with the kernel K. It oers new features, in particular the computation of the maximal ball centered at a point of the polygon. Such a ball has a center which is a Point_3 from K, and a squared radius whose type is NT.

The candidate balls. Following the results of section 4.1, the candidate balls used are only centered on the vertices of the medial axis. Such balls are associated with the medial axis, as just discussed.

The volume of the selected balls. Computing the volume of a union of balls is a dicult problem, from a combinatorial, but also numerical standpointinverse trigonometric functions are involved. We use our certied algorithm [START_REF] Cazals | Computing the volume of union of balls: a certied algorithm[END_REF] which returns an interval certied to contain the exact volume. More precisely, due to the impossibility to obtain a volume as an exact number type, whatever the kernel used (exact, inexact), the centers and radii of the candidate balls are converted to doubles. These balls are input to our algorithm, which requires two template parameters: the number type of the output (double or interval), and the level of exactness used to compute the constructions involved in the volume computation, namely the coordinates of Voronoi vertices, and boundary points of the union of the selected balls. Following the discussion in [START_REF] Cazals | Computing the volume of union of balls: a certied algorithm[END_REF], the three options are referred to as (faster, ck_pt_exact and all_exact). Practically, we use the pair (double, faster) for the inexact kernel, and (interval, all_exact) for the exact kernel.

Results

Dataset.

As test set, we used the 96 protein -protein complexes from [START_REF] Lo Conte | The atomic structure of protein-protein recognition sites[END_REF]. The complexes are of high biological interest (all of them are coupled to well identied biological processes). The number of atoms lies in the range [1008, 13214], with a median of 3757.

Performances and robustness issues.

The properties of predicates and constructions of the exact and inexact kernels has been discussed in section section 5.2. We compared the volume ratios obtained with these two options on a set of 10 protein complexes, and did not observe any dierence before the third digit. For two examples discussed in detail below, these ratios are ∼ 0.69 for the protein complex 3sgb with r = 2.8 and 20 selected balls, and ∼ 0.64 for the immunoglobulin 1igt with r = 2.8 and 103 selected balls. For running times, we compared the execution time for the construction of DT B, for DT V and for the medial axis. The selection itself was excluded since from the timing, as also noticed in section 5.2, our volume computation algorithm uses double as number type. On the aforementioned 10 models, we observed that the exact kernel was on average about 150 times slower than the inexact one. For these two reasons absence of obvious degeneracies and much better running time, the results reported in the sequel were computed with the inexact kernel.

Using the inexact kernel, it is observed that the running times for computing DT B and DT V are a mere order of magnitude slower than the CGAL ones3 for the regular triangulation case (Fig. 5). These running times are naturally consistent with the fact that the geometric objects manipulated behave nicely for our molecular models: both the number of boundary points (Fig. 6) and the primitives of the medial axis (Fig. 7) are linear in the number of input balls.

Inner approximation guarantees.

Intuitively, the ability of greedy to provide a good approximation relies on the possibility to choose large balls, which depends itself on two parameters. First, the topological complexity: the closest to a topological ball the domain F C , the better. Second, the geometric complexity: the more convex the domain F C , the better. Before commenting these properties on molecular systems, recall that in a vdW model, the radii of the balls vary in the range 1-2Å, and that only balls of atoms linked by covalent bonds intersect. Thus, for a vdW model, one expects the volume covered to vary linearly as a function of the selection size, which is exactly observed (Fig. 8, red curves). Now, enlarging the input balls by a quantity r w , e.g. r w = 1.4 to dene a solvent accessible model, results in simplifying the topology of F C . Thus, the larger r w , the larger the candidate balls, and the better the ratio curve (Fig. 8

again).

For a xed budget of balls and a given expansion radius, e.g. r w = 1.4, one also expect complex topologies, characterized by a high Euler characteristic, to yield more dicult problems. This trend is also observed (Fig. 9).

As for the incidence of the overall shape, a roughly convex system (Fig. 10(A) versus Fig. 11(A)) clearly yields more favorable volume ratio curves (Fig. 8(Top, Bottom)).

Coarse graining molecular models.

The immunoglobulin (Ig) structure just used provides a good example, because of its non convexity, to test the algorithm to go beyond the inner covering problem. This model also exhibits a topological rather than geometric diculty, since at the tip of the two arms each so-called variable domain has the topology of a lled torus (Fig. 11(A)). To this end, in a manner identical to the scale axis transform [START_REF] Giesen | The scale axis transform[END_REF], we may enlarge the model by r w , approximate it, and compare the balls obtained with the initial vdW model. When implemented with r w = 2.8, this heuristic yields a coarse-grain model exhibiting a discrepancy of about 3 atoms sticking out from the coarse grain model in the worst-case (Fig. 11(C)). Moreover, the topology of the variable domains has also been preserved. Thus, one gets an accurate model for a number of balls which has been divided by about 100 in this case. (In biophysics, a resolution of 5Å is the maximum which may be tolerated to say that a model has atomic resolution.)

Conclusion and Outlook

This paper studies two basic problems dealing with collection of balls, namely selecting a diverse set of molecular conformations, and providing an accurate inner approximation of a molecular model. Both problems are shown to be geometric versions of max k-cover, the weight function being a function of the geometry of the cells of a surface or volumetric arrangement, rather than being uniform as in the combinatorial setting. Yet, for the volumetric case, the approximation bound known in the combinatorial setting is preserved, provided that the weights are nonnegative. The implementation of our algorithms hinges upon state-of-the-art software coupled to the CGAL library. In particular, this implementation involves the exact calculation of a Delaunay triangulation for points whose coordinates are degree two algebraic numbers, the intersection of the dual of this triangulation with the α-complex of the input balls, and the certied calculation of the volume of a union of medial balls. This implementation handles molecular models containing up to O(10 5 ) atoms within minutes. For these reasons, we believe that our algorithms, the inner approximation in particular, should prove useful for a broad class of geometric approximation problems dealing with balls, in particular in the context of approximate medial axis transforms, where the focus has been so far on asymptotic propertiesupon increasing the number of balls. Yet, our work calls for further developments, both in the theoretical and applied directions. On the theoretical side, two challenging questions are of high interest. First, our greedy algorithms come with guarantees for the inner approximation problem. But monitoring the symmetric dierence of the input domain and of the selection, or the Hausdor distance between their boundaries is also clearly of interest. Second, constraining the geometric selection by topological criteria, e.g. prescribed Betti numbers, would also be of the highest interest. However, approximation problems aiming at accommodating both geometric and topological criteria are likely to be challengingit has been shown that the so-called homology localization problem is NP-hard. In an applied vein and as mentioned in introduction, we believe that a key application of our algorithms will be the design of coarse-grain macro-molecular models, to investigate macro-molecular machines and simulate crowded environments within whole cells. But prior to undertaking these challenges, one will have to decorate our purely geometric coarse-grain models with bio-physical properties.

RR n°8205

Figure 3 Greedy selection may not yield the optimal solution. Greedily selecting two sets out of C 1 , . . . , C 5 yields a score of 12 (selecting C 3 and then C 2 ), while the optimum is 14 (selecting C 4 and C 5 ). The shaded cells have the weights as indicated and the unshaded cells have null weights.

A 5 4

A 6 4 A 3 2 A 4 2 A 1 1 A 2 1 C 3 8 C 2 4 C 1 2 C 4 7 C 5 7

  4.1, one has B * d ⊂ B * e ∪ B * f and similarly B * p ⊂ B * c ∪ B * d . Thus, B * p ⊂ B * c ∪ B * e ∪ B * f .

Figure 1 Figure 2

 12 Figure 1 2D conformations, each consisting of four ballsrst and fourth balls are common, and the induced surface and volume arrangements. The (two dimensional) volume occupied by the two conformations is decomposed into 19 cells (boxed numerals). The circled numerals feature the surface arrangement of the ball centered at a 1 , based on intersections with neighboring balls.

Figure 4 AFigure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10

 45678910 Figure 4 A tight example for the greedy strategy.

Figure 11

 11 Figure 11 Coarse graining a molecular model: the example of an immunoglobulin (PDB id: 1igt). (A) The atomic van der Walls models. (B) A coarse grain model of 103 balls, dened as the inner approximation of the atoms whose van der Waals radii have been expanded by r w = 2.8. Three balls were added to the greedy selection with s = 100 as explained in section 5.1 to force the connectivity. (C) The superimposition of both models. Note in particular that the hole of the arm on the left hand side is respected. 12,543 atoms 103 balls (A) (B)

RR n°8205

http://mgl.scripps.edu/people/goodsell/books/MoL2-preview.html RR n°8205

http://www.cgal.org Inria

http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Triangulation_3/Chapter_main. html#Subsection_39.6.1 Inria

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments. Michael Hemmer is acknowledged for his help with the CGAL kernels. This work has partially been supported by the Computational Geometric Learning STREP project of the EC 7th Framework Programme (EC contract No. 255827).


Gilles Rabatel
email: gilles.rabatel@irstea.fr

Nathalie Gorretta

Sylvain Labbé

Getting simultaneous red and near infrared bands from a single digital camera for plant monitoring applications

Keywords: multispectral sensor, near-infrared band, NDVI

Multispectral images including red and near-infrared bands have proved their efficiency for vegetation-soil discrimination and agricultural monitoring in remote sensing applications. But they remain rarely used in ground and UAV imagery, due to a limited availability of adequate 2D imaging devices. In this paper, we propose and evaluate an original solution to obtain simultaneously the near-infrared and red bands from a standard RGB camera, after having removed the near-infrared blocking filter inside. First, the theoretical approach is described, as well as simulated results on a set of soil and vegetation luminance spectra with two different still cameras (Canon 500D and Sigma SD14). Then examples of images obtained in real field conditions are given, and compared with standard colour image acquisition for pixel-based plant/soil discrimination, using an automatic thresholding method. It appears that in most cases our new acquisition procedure brings a significative improvement, opening new opportunities for crop monitoring applications.

Introduction

The Normalized Difference Vegetation Index, or NDVI, introduced in the early seventies by [START_REF] Rouse | Monitoring vegetation systems in the great plains with ERTS[END_REF], remains today a very popular tool in the remote sensing community dealing with agricultural monitoring. This is mainly due to its remarkable ability to discriminate vegetation from other material in multispectral satellite images.

Green vegetation is characterized by a high reflectance in the near-infrared domain (typically 50 to 80%), which contrasts with a very low reflectance in the red wavelengths, due to chlorophyll absorption. Let us call R and NIR the digital counts obtained through the red and the near infrared bands of a multispectral sensor. The NDVI, expressed as: NDVI = (NIR-R)/(NIR+R)

(1) is a scalar value in the range [-1, 1]. The higher is this value, the higher is the probability that it corresponds to green vegetation. By extension, numerous attempts have be made to directly link NDVI or other derived indexes based on R and NIR to agronomical indices such as biomass and LAI or Leaf Area Index [START_REF] Jindong | Assessing broadband vegetation indices and QuickBird data in estimating leaf area index of corn and potato canopies[END_REF].

The popularity of NDVI in remote sensing has been widely supported by the availability of R and NIR channels on most satellite line-scanning sensors (Landsat, SPOT, etc.). Unfortunately, it is not the case for crop monitoring applications at a lower spatial scale: most vision systems embedded on ground vehicles or UAV (Unmanned Aerial Vehicle), which require 2D sensors, are based on standard colour cameras, leading to robustness issues in vegetation detection.

To face this situation, some camera end users requiring a near infrared channel have developed alternative solutions around standard RGB cameras, taking benefit of an "undesirable" property of their silicium sensing array: because the colour filters in the Bayer matrix have a filtering action limited to the visible domain, the camera manufacturers are constrained to add a near-infrared blocking filter to match natural colorimetry requirements.

By removing this additional filter, we obtain a "modified" camera sensitive to near infrared wavelengths. A first possibility to get separate R and NIR bands is thus to use simultaneously a standard and a modified colour camera, the second one being equipped with a near infrared pass-band filter. However, important issues arise concerning the pixel alignment of the two images obtained [START_REF] Lebourgeois | Can Commercial Digital Cameras Be Used as Multispectral Sensors? A Crop Monitoring Test[END_REF].

Another interesting approach is to use a single modified camera, and to recover the red and near infrared bands as specific linear combinations of the three raw channels, provided an adequate filter is settled in front of the camera lens. The theoretical base of this approach has been detailed by [START_REF] Rabatel | Getting NDVI Spectral Bands from a Single Standard RGB Digital Camera: A Methodological Approach[END_REF]. In the present paper, we will just recall its main results, and then emphasize on the results obtained in real field conditions.

For this purpose, various aerial images of field crop have been acquired, using both modified and standard still cameras, and their ability to discriminate between soil and crop have been compared.

Material and methods

Red and near infrared recovering from a single camera

Let us consider a three-channels camera, where the channels are characterized by their respective spectral sensitivities c 1 (O), c 2 (O), c 3 (O), and a virtual channel sensitivity s(O) expressed as their linear combination :

s(O) = a 1 .c 1 (O) + a 2 .c 2 (O) + a 3 .c 3 (O)
. s(O) can be easily simulated from the camera, the digital count for a given pixel being equal to:

S = a 1 .C 1 + a 2 .C 2 + a 3 .C 3 (2)
where C 1 , C 2 , C 3 are the raw digital counts for this pixel In that sense, c 1 (O), c 2 (O), c 3 (O) build a 3D subspace of the general spectral space (with infinite dimension). Every possible simulated sensitivities are contained in this 3D subspace. Now let us consider a desired virtual sensitivity v(O) (i.e. near infrared band). In the general case, v(O) will not be included in the subspace generated by (c 1 (O), c 2 (O), c 3 (O)), and the best approximation will be its projection P(v) on this subspace (Fig. 1).

(c 1 , c 2 , c 3 ) subspace v P(v) The quality of the approximation can be evaluated by the angle between the vectors v and P(v), known in spectrometry as Spectral Angle Mapper or SAM [START_REF] Yuhas | Discrimination among semi-arid landscape endmembers using the spectral angle mapper (SAM) algorithm[END_REF].

Obviously, a standard RGB camera will allow recovering properly the red band (it has been designed for it) but not a near-infrared one, because the corresponding wavelengths (or dimensions) have been removed. What about a modified camera with no blocking filter ?

The approach involved here consists in searching for the best possible approximation for both red and near infrared bands, using an additional filter in order to modify the c 1 (O), c 2 (O), c 3 (O) sensitivities, and thus the 3D projection sub-space. For technical reasons, only lowpass filters (cutting the short wavelengths) have been considered.

This study has been led with two high-level SLR (Single lens Reflex) cameras: Canon 500D (Canon, Tokyo, Japan) and Sigma SD14 (Sigma, Kawasaki, Japan), this last one having the particularity of being built around a 3-layered CCD sensor, instead of using a traditional matrix of microfilters. For both cameras, the near infrared blocking filter has been removed, and the resulting sensitivity curves c 1 (O), c 2 (O), c 3 (O) have been measured in laboratory using the tunable monochromatic light source of a laboratory spectrometer (V-570, Jasco Inc, Easton, USA). Projections of standard Landsat red and near infrared bands have then been simulated for various low-pass filter cutting wavelengths using Matlab 7 (The MathWorks, Natick, MA, USA). For each cutting wavelength, the quality of the projections Q has been computed as the sum of SAM (in radian units):

Q = SAM R + SAM NIR = SAM(r(O), P(r(O))) + SAM(nir(O), P(nir(O))) (3)
where r(O) and nir(O) are respectively the targeted red and near infrared sensitivities.

Finally, NDVI values have been computed for a set of vegetation and soil luminance spectra, and compared with NDVI issued from standard red and near infrared bands. These spectra were issued from hyperspectral images of wheat durum acquired in experimental fields (INRA, Domaine de Melgueil, France) in march 2011, using a Hyspex VNIR-1600 camera (Norsk Elektro Optikk, Norway).

Field crop image acquisition

According to results of the study above, various aerial images have been acquired using modified cameras (no infrared blocking filter) equipped with a Rot R8 (25 A) glass filter (Hama GmbH & Co, Monheim, Germany) in front of the lens. Depending on their availability, either a Canon 500D or a Canon 350D was used. Images were acquired by a private company (AvionJaune, Montferrier, France) using planes or Unmanned Aerial Vehicles (UAV), depending on the flight altitude.

For comparison purpose, the same aerial scenes have also been acquired using standard RGB cameras. In order to illustrate the plant/soil discrimination ability in both cases, the following procedure has been followed:

-for images acquired with the modified cameras, NDVI values have been computed using the R and NIR digital counts obtained from the procedure described in 2.1.

-for standard RGB images, the Excessive Green Index EGI = 2G-R-B as well as its normalized version NEG = EGI/(R+V+B) [START_REF] Woebbecke | Color indices for weed identification under various soil, residue, and lighting conditions[END_REF] have been computed.

-NDVI, EGI and NEG images have then been computed by scaling the corresponding index values in the range 0-255, and thresholded according to their grey level histogram. For this purpose, we used the Otsu method [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF], which minimises the ratio between inter-class inertia and total inertia for the two classes issued from thresholding. This ratio, known as Wilk's lambda (/ w ) in multivariate analysis, gives an indication on the thresholding quality and robustness.

In the following examples, NEG results will not be presented, because their / w values were always lower than those obtained with EGI.

Results and discussion

3.1. Red and near infrared band projections Fig. 2 shows the quality Q, as expressed in equation ( 3), of red and near infrared band projections for Canon 500D and Sigma SD14 cameras when varying the long-pass filter cutting wavelength. FIGURE 2: Red and near-infrared projection quality (in radians) as a function of the cutting wavelength It appears that in both cases, the optimal projection is obtained for a cutting wavelength of about 600 nm, which corresponds to a red low pass filter. However, the best result is obtained with the Canon camera. For this reason, further studies will be limited to this type of camera, equipped with a photographic filter Hama Rot R8 (25 A), cutting wavelengths below 600 nm. The corresponding red and near-infrared projections are expressed as:

P(r(O)) = 0.0968* c 1 (O) -0.1722* c 2 (O) + 0.0842* c 3 (O) (4) P(nir(O)) = -0.0455* c 1 (O) + 0.0010* c 2 (O) + 0.2605* c 3 (O) (5)
Fig. 3 shows the shape of these projections, as well as the corresponding NDVI values computed on our set of plant and soil spectra. The better quality of NDVI thresholding compared to EGI thresholding is confirmed by the Wilks lambda values (resp 0.80 and 0.70). Also, we can notice that equations ( 5) and ( 6 According to images and Wilks lambda values, the EG index gives a better thresholding result. However, only the trees are discriminated from the soil. The vine rows (dark green) are not detected. In that sense, NDVI is more robust for plant detection.

Conclusion

We have shown in this paper that NIR and R bands can be issued from a single standard digital RGB still camera, by replacing the internal near-infrared blocking filter by an off-theshelf low-pass filter set in front of the lens. Though the resulting red and near-infrared sensitivities do not match exactly the standard curves r(O) and nir(O), they provide satisfactory NDVI values for plant and soil, and thus bring more robustness in plant-soil discrimination, compared to standard color images.

These results are particularly interesting for in-field or low altitude imaging applications, which are presently mainly limited to standard color acquisition devices. In that sense, they open new possibilities in terms of high-resolution imaging for crop monitoring.

FIGURE 1 :

 1 FIGURE 1: Approximation of a desired spectral sensitivity by projection

FIGURE 3 :FIGURE 4 :

 34 FIGURE 3: left: red and near-infrared projections P(r(O)) and P(nir(O)) using the Rot R8 filter (in dot lines: targeted sensitivities r(O) and nir(O)). right: NDVI values computed with r(O), nir(O) (abscissa) and P(r(O)),P(nir(O)) (ordinate) on a set of soil and plant spectra

  FIGURE 5: EGI and NDVI computation on a vine parcel a,b: color image (full scale and detail); c,d: EGI computing and thresholding (/ w =0.7333); e,f : modified camera NDVI computing and thresholding (/ w = 0.70)

CIGR-Ageng 2012. International Conference on Agricultural Engineering, 8-12 Juillet 2012, Valence, Espagne

Acknowledgements

The research leading to these results has received funding from the European Union's Seventh Framework Programme [FP7/2007[FP7/ -2013] ] under grant agreement n°245986.

We acknowledge L'avion jaune (http://www.lavionjaune.fr/) for having provided the modified cameras and made the aerial image acquisitions.


published or not. The documents may come

The complexity in the French Wood Supply Chain reduction: A traceability perspective

T. Saikouk,


Jean-Rémy Chardonnet

Jean-Claude Léon

Designing Interaction in Virtual Worlds through a Passive Haptic Peripheral

This paper presents a prototype of a hands-on immersive peripheral device for controlling a virtual hand with high dexterity. Based on the results of users' tests on previous versions of our device and on the analysis of a manipulation task, this prototype is as easy as a mouse to use and allows the control of a high number of degrees of freedom (dofs) with tactile feedback. Design issues, physical phenomena and physiological behaviors are tightly linked and highly influence interaction. The goals corresponding to these issues include the choice of sensors' technology and their position on the device, low efforts exerted while using the device, relevant multisensorial feedback, performance of achieved tasks. An example of a grasping task illustrates the effectiveness of our device to achieve intuitive and efficient interactions, bringing new insights for collaborative interaction.

I. INTRODUCTION

Object manipulation is a natural and essential everydaylife task that can be done alone or in collaboration with other persons or robots. Manipulating an object is conditioned by the number and the quality of information returned to the brain. This paper presents an extension of the HandNavigator, an immersive peripheral for controlling a virtual hand in virtual environments [START_REF] Chardonnet | Design of an immersive peripheral for object grasping[END_REF] for precise object manipulation.

A. Related work

Interacting with an object is a difficult task to achieve since several possible grasping configurations are possible. Thus, imitating a grasping motion in virtual environments is not easy. This diversity of configurations, conditioned by several parameters such as the shape, is not addressed by the current devices and software. Indeed, each movement of a hand and its fingers cannot be controlled in virtual environments because the interfaces complexity dramatically increases. Nevertheless, some natural configurations of a hand can be achieved as long as interfaces give enough immersive sensations. Through an appropriate peripheral device allowing some motions of a hand, it is possible to evolve from a qualitative validation of a task to a quantitative one where hand configurations can be formalized and quantified.

Several solutions have been proposed to perform general manipulation tasks in virtual environments, transforming motions of a real hand to motions of a virtual one [START_REF] Bowman | 3D User Interfaces: Theory and Pratice[END_REF]. Often, corresponding solutions rely on motion capture based systems, for instance optical systems through cameras [START_REF] Dewaele | Hand motion from 3d point trajectories and a smooth surface model[END_REF], 1 J.-R. Chardonnet is with Arts et Métiers ParisTech, CNRS, Le2i, Institut Image, 71100 Chalon-sur-Saône, France jean-remy.chardonnet at ensam.eu oratory, 38041 Grenoble Cedex 9, France jean-claude.leon at grenoble-inp.fr [4], or mechanical systems through data gloves where sensors are attached to an exoskeleton covering the hand [START_REF] Sturman | Hands-on interaction with virtual environments[END_REF]. These solutions allow a high number of possible configurations of a hand, thus generating multiple possibilities to move a hand at real scale around the objects and, consequently, to grasp and manipulate them. However, they have several major drawbacks. These systems require preliminary calibration, which can be time consuming. Their capability to generate a large range of motions drastically increases their mechanical complexity and hence, their cost. As mentioned above: calibration is complex but using such systems is also difficult. Consequently, these systems are not ready-to-use. More precisely, it is difficult to load pre-saved configurations or to analyze them for issues related to accessibilities or trajectories, e.g., for collaborative tasks. Also, it is not possible to interrupt, at any time, a current task without losing its parameters and current configuration. Contact is hard to interpret even with some feedback, especially in the case of continuous contact, e.g., sliding a finger on a surface without loosing contact. Optical based systems suffer from occlusion problems. Indeed, there is often a point that is not visible by the camera. To solve this issue, cameras can be added but cost and space highly increase without avoiding configurations where occlusion problems remain. From an ergonomic point of view, achieving long tasks with a raised hand configuration quickly generates a muscular fatigue.

Simple systems, designed to be integrated in desktop environments, have been proposed, such as the SpaceNavigator, to achieve 3D interactions. They allow a maximum of six degrees of freedom (dofs) and a comfortable use because the user does not use raised hand configurations. These devices are well suited for pointing and object moving tasks. Their main limitation comes from their number of dofs that is not appropriate for more complex manipulation tasks requiring more dexterity, as required for example for col-laborative tasks with other avatars (humans, robots). Multitouch pads are more and more used in laptop computers, smartphones and allow more complex interaction with virtual environments [START_REF] Rosenberg | The unmousepad -an interpolating multitouch force-sensing input pad[END_REF], [START_REF] Hook | A reconfigurable ferromagnetic input device[END_REF]. Their main advantage lies in their high sensitivity, allowing a good ergonomics for long tasks (from tens of minutes to one hour). However, they are designed mainly for 2D applications. Recently, some 3D applications are emerging but they cannot be enough immersive because these interactions are still performed in 2D.

Some systems integrate a haptic feedback [START_REF] Kim | Development of tension based haptic interface and possibility of its application to virtual reality[END_REF], [START_REF] Bouzit | The rutgers master ii-new design force-feedback glove[END_REF]. This feedback can also be found on larger equipment, such as haptic devices used for instance to interact in physical simulations (see for example [START_REF] Allard | Sofa -an open source framework for medical simulation[END_REF], [START_REF] Chardonnet | Interactive dynamic simulator for multibody systems[END_REF]), and includes a force feedback. These solutions increase the amount of information sent to the user and thus allow to better take advantage of the sensory capabilities of his/her brain. Haptic feedback is perceived only when the user touches a virtual object but is lacking totally or partially when moving freely in 3D space, without any contact [START_REF] Johansson | Sensory input and control of grip[END_REF], [START_REF] Zhai | The influence of muscle groups on performance of multiple degree-of-freedom input[END_REF]. Active haptic systems' complexity increases as the number of force feedback components increases. These systems allow a user to better feel the virtual objects, but are hardly used by the general audience as they are expensive and they have a technological complexity that needs a good knowledge of the device. Moreover, one major limitation is parameters tuning, required to get realistic sensations of physical models, which is time-consuming and purely subjective.

One last solution is to provide passive feedback. This feedback can be perceived as a improvement compared to systems without any feedback [START_REF] Insko | Passive haptics significantly enhances virtual environments[END_REF], and even if it does not allow a user to get a total immersion in virtual environments, it can fool the user's proprioceptive senses [START_REF] Lecuyer | Pseudo-haptic feedback: can isometric input devices simulate force feedback?[END_REF] at low cost by using simple objects such as sponges or vibrators [START_REF] Zhai | The influence of muscle groups on performance of multiple degree-of-freedom input[END_REF], [START_REF] Pai | The tango: A tangible tangoreceptive whole-hand human interface[END_REF]. It can easily be integrated in classical peripheral devices.

B. General presentation

Regarding the approaches above, our solution is based on a passive haptic feedback and a hands-on interaction. We present an extension of the HandNavigator described in [START_REF] Chardonnet | Design of an immersive peripheral for object grasping[END_REF]. Especially, validation tests of the existing prototypes are conducted by several users on grasping tasks scenarii, leading to the design of a new requirements of the HandNavigator with a more ergonomic shape, sensors allowing better dexterity and interactions, while integrating tactile feedback for an enhanced immersion, which is not the case in [START_REF] Chardonnet | Design of an immersive peripheral for object grasping[END_REF].

In this paper, we will at first expose the issues related to different possible configurations for object manipulation. Then, we introduce the previous prototypes and their results of validation tests performed on these prototypes, these results are analysed for the design of a new version. The following section describes the new prototype and shows a simple example using this prototype before concluding.

II. ISSUES AND GOALS

The kinematic structure of our hand allows us a wide range of configurations, among which we select the best one to perform the desired task. This ability to adapt the configurations to objects greatly improves our dexterity. Dexterity reduces here as the coordination between a hand, its fingers and the user's eyes, implying the ability to use them to perform activities requiring high accuracy. Thus, visual and tactile feedbacks are mandatory to achieve fine dexterity [START_REF] Talati | Visual and tactile guidance of dexterous manipulation tasks: an fmri study[END_REF]. More precisely, the concept of dexterity here is bounded by a group of tasks we want to perform with our peripheral device. These tasks can be summarized as follows:

• contact-free motion of each virtual finger independently or not from the others, as natural as possible;

• manipulate objects as naturally as possible. To achieve them, we must focus on the kinematics allowing a user to control independently the global motion of his virtual hand and that of its fingers. The simultaneous mobility of a hand and its fingers is mandatory to perform a large panel of tasks when a user interacts with the device. When analyzing a real grasping motion, we can distinguish three steps:

1) the hand and its fingers move freely in 3D space to reach the object to be grabbed. During this step, the user already adapts the configuration of his hand to fit the object's shape; 2) the hand and its fingers touch the object. Contact can appear sequentially (the hand moves towards the object then its fingers bend to touch it) or progressively (the hand moves while closing its fingers). During this step, tactile and visual feedbacks are necessary and guide the user while performing this task; 3) the user applies a force on the object to grip it. This force depends on the friction between the object and the hand and its fingers. For instance, the low friction of slippy objects requires higher gripping forces so that the objects do not slip. Here also, kinesthetic and touch feedbacks guide the user to act on the object in the best way. During this step, his hand configuration evolves, thus generating different hand's postures. To perform these steps, the user controls a high number of dofs. In virtual worlds, several dofs of a virtual hand mechanism can be configured using inverse kinematics, allowing a user to generate a large number of postures while reducing the number of dofs to control on the device. However, it is mandatory to monitor more than six dofs (defining the global hand position and orientation) to perform precise tasks with high dexterity and one dof per finger is a minimum.

Considering these characteristics, the device we develop must take into account the following requirements: (i) allow the user to control simultaneously a high number of dofs, (ii) allow a user to perform complex hand and finger motions, (iii) ergonomics to avoid pain or fatigue that drastically reduces motion dexterity and does not allow a long use (from several minutes to several hours), (iv) return to the user relevant contact information, and (v) be cheap, easy and ready to use, and calibration-free to ease its integration.

III. ANALYSIS OF THE EXISTING PROTOTYPES A. Existing prototypes and users' tests

We describe the main features of different prototypes of the HandNavigator and users' tests performed on these proto- types. The HandNavigator consists in two major components: a 3D mouse SpaceNavigator 1 and, on top, a housing where sensors for fingers are located. The prototypes presented in [START_REF] Chardonnet | Design of an immersive peripheral for object grasping[END_REF] use either pressure sensors, trackballs, or scrollpads (see Fig. 2). We refer the reader to [START_REF] Chardonnet | Design of an immersive peripheral for object grasping[END_REF] for a detailed description of the SpaceNavigator and the different versions of the HandNavigator. We call these prototypes V2, V3t, and V3s respectively, according to their chronology.

To measure the performance, the usability, the controllability of each prototype and validate issues linked to human perception and learning capabilities, we performed tests with unexperimented users on each prototype. With each prototype, each user followed the procedure hereunder:

1) the user manipulates the devices for three minutes. Two sub-tasks are proposed:

• firstly, he concentrates on the motion of the virtual hand and tries to reach a green sphere fixed in the virtual scene. A yellow sphere, attached to the virtual hand, helps the user visualizing a target area on this object. When the two spheres meet each other, the green one changes color to indicate that the global hand position is satisfactory;

• then, the user concentrates on the motion of the virtual fingers. The virtual hand stands still in the virtual scene to allow the user acting freely on the sensors to monitor the virtual fingers, without generating any perturbation on the SpaceNavigator in charge of the global hand position; 2) subsequently, the user grasps a virtual giraffe located at some position in the virtual scene (see Fig. 3). The user gets feedback when the giraffe is grasped, since it changes color. The user must repeat this task several times. During each trial, the time required to achieve this task is measured to get a learning curve of the user's adaptation to the different prototypes. At the end of the test, users had to fill up the questionnaire reproduced in Table I. The users score each criterion from 0 to 4. Note that, compared to other questionnaires of the literature, e.g., [START_REF] Wiese | Investigating the learnability of immersive free-hand sketching[END_REF] or the NASA-TLX protocol [START_REF] Hart | Development of nasa-tlx (task load index): Results of empirical and theoretical research[END_REF], devoted to task achievement evaluation, we adapted their content to evaluate our prototypes. We asked 32 users aged between 20 to 25, right-handed, left-handed, men and women, to perform these tests. We obtained the results depicted in Fig. 4.

We observe that the prototypes get better scores as they evolve, especially in terms of comfort and usability, the users feel less pain or fatigue with the last prototype V3s, with 1 http://www.3dconnexion.com 4)) 2. Motion of the fingers (difficult (0) -easy (4))

3. Motion of the fingers without moving the virtual hand (difficult (0) -easy (4)) [START_REF] Schlattman | Simultaneous 4 gestures 6 dof realtime two-hand tracking without any markers[END_REF]. Comfort (bad (0) -good (4)) 5. Finger positions with respect to the sensors (bad (0) -good (4)) 6. Sensors' monitoring (difficult (0) to easy (4)) 7. Size of the prototype (inadequate (0) -well fitted (4)) improved performance. However, we found out that the users had difficulties to move the virtual fingers without moving the virtual hand because they held the device too tightly. We think that it is a psychological phenomenon due to the novelty of the device that causes stress, exactly in the same way people using a computer mouse for the first time. An extended use shows that this stress tends to disappear.

The learning curve for each prototype is depicted in Fig. 5. We observe an adequacy between the qualitative results in Fig. 4 and the required completion time. The prototype V3s is easier to control because after only three trials the users can grasp the giraffe in less that 5s whereas with the first prototype: V2, the time does not reach less than nearly 20s. Note that the prototype V3s seems to be the most intuitive and the one the users are immediately the most familiar with.

B. Discussion and new design requirements

The users' tests teach us several important issues allowing us to design a new prototype that is more efficient in terms of interactivity, dexterity and comfort.

The first issue is technological and addresses the sensors used to control the finger motion. To be efficient and avoid user's fatigue and undesired motions, we should consider a range of sensors allowing continuous motion of the real fingers between the open and closed configurations of the virtual fingers. One of the main reasons for which prototype V2 is badly evaluated by the users is the latency in finger motion: there are two sensors for each finger, one to open and one to close it. V3t uses trackballs and does not exhibit latency when opening/closing the real fingers because the real fingers just rotate slightly the sensor to switch from virtual finger's opening to closing. However, the user must roll them several times to reach the desired finger configuration, which, for long tasks, can be painful. V3s uses scrollpads and allows the user moving continuously and smoothly the virtual fingers. The major difficulty with these sensors is their constraint set on the user's fingertips having to travel a long distance to generate large virtual finger motions, increasing the probability of producing undesired motions. As a conclusion, it appears important to identify sensors capable of producing a high ratio of virtual finger to real finger displacement. Therefore, we will consider another sensor's technology: lever-switches, that can be found for example on digital cameras to zoom in or out. If the user releases the sensor, it goes back to its neutral position and the virtual finger stops.

The second issue relates to neuro-psycho-physiology. The interest of relying on small motions of real fingers to produce large motions of the virtual fingers is not only mechanical: it avoids solliciting too much the user's muscles of his arm and forearm. We observed during the evaluation tests that users were stressed on the prototypes. When reducing their muscular activity, we ensure them a better dexterity and performance, coupled with better comfort of use. Besides, generating full scale motions of the hand and the fingers is useless, as obtained with data gloves for example. On the one hand, with the sensors' technologies used, we can get a fine enough dexterity and accuracy. On the other hand, several neuro-psychological studies show that for a same task, between motor activity (the effective task) and mental activity (the pure cerebral activity of a task), very small difference is observed in terms of duration [START_REF] Kosslyn | Mental rotation of objects versus hands neural mechanisms revealed by positron emission tomography[END_REF], [START_REF] Harris | Selective right parietal lobe activation during mental rotation -a parametric pet study[END_REF]. In other words, if the user, without moving his hand and fingers, performs a pure cerebral task, he can get roughly the same performance results as if he is really performing the motion while producing a very small muscle activity. Finally, other studies show that for tasks requiring high precision and dexterity, the hand is the most efficient limb of the human body. Indeed, the fingers, more specifically the thumb, and the hand occupy the largest area of the somatosensorial cortex compared to other limbs [START_REF] Penfield | The Cerebral Cortex of Man: A Clinical Study of Localization of Function[END_REF], [START_REF] Balakrishnan | Performance differences in the fingers, wrist, and forearm in computer input control[END_REF]. Thus, it is important to focus specifically on the way to use the device to increase performance and dexterity.

The last issue focuses on the tactile feedback when contact occurs between the virtual hand and an object. Indeed, just looking at an object is not enough, especially in virtual worlds, to determine the way to grasp it because parameters such as friction, temperature cannot be sensed if there is no contact, and occlusions can appear. In the evaluation tests, none of the prototypes include a tactile feedback, which can explain the time required to grasp the giraffe (at least 5s) whereas it is initially close (20cm) to the virtual hand. In real conditions and at equal displacement speed, it should take only one second to grasp it. Therefore we will set up a passive haptic feedback, enough to fool the users's proprioceptive senses.

We see that the design of an immersive device goes beyond pure mechatronics and must also consider more general issues linked to neuro-psycho-physiological activities.

IV. NEW PROTOTYPE STRUCTURE

A. Individual module

In order to incorporate the different requirements mentioned above, we designed a module allowing the user to control just one finger. This module is depicted in Fig. 6. It contains a lever-switch, for controlling a virtual finger with a low motion amplitude of a real finger, a pressure sensor and a vibrator to give the user touch and grasp feelings. The lever-switch is active only when the finger is free to move in 3D space, i.e., without any contact. The virtual finger is velocity-based controlled, allowing a user to interrupt and resume a task anytime without loosing its configuration. The vibrator is activated only when the virtual finger touches a virtual object and vibrates for a short time (half a second). The pressure sensor aims at giving the sensation to the user that he really grasps a virtual object and it is active only while the user is holding it. Indeed, the pressure sensor has a lower sensitivity than the other sensors we incorporated [START_REF] Chardonnet | Design of an immersive peripheral for object grasping[END_REF], as the user must apply higher forces (1N) than to the SpaceNavigator (0.4N) to activate it. Through this intrinsic characteristic, the user has to apply a significant force, similar to a configuration he encounters when holding tight a real object. Finally, a spring between the lever-switch and the pressure sensor separates the two modes (free and constrained hand motions). Finally, since one of our requirement devotes the device to a general audience, each module is inserted in an armhole so that the module can be adapted to a large range of hand sizes. The interest of separating the module from the housing is not only ergonomics, it enables setting up vibration barriers to avoid their propagation in the whole structure, which will be of major concern for the final prototype.

B. Complete prototype

We duplicate the module to include the other fingers. We completely redesigned the shape of the prototype to be ergonomic and allow better interaction. Similarly to the previous prototypes and for the same reasons related to the hand kinematics, it is possible to control only four fingers. Each finger has one sliding module for different hand sizes. When a user touches an object, he can immediately know which finger touches the object, thanks to individual vibrators. The final prototype, called V4, is depicted in Fig. 7. The schematic design of the prototype is shown in Fig. 8. As for the previous prototypes, the SpaceNavigator allows the control of the 6 dofs of the virtual hand's wrist. Each module controls the joint position of the last phalanx of a virtual finger. The joint position of the other phalanxes are equal to the one of the last phalanx to simplify the device. We used laser sintering techniques to manufacture this prototype and added a silicon-based coating as damping material to filter out the vibrations and isolate each finger.

V. DEMONSTRATION

We now show an evaluation of prototype V4 through a task example. The goal is to position the virtual hand and fingers on a part coming from a food processor as a manual assembly task. As shown in Fig. 9, there are several possible configurations to grasp it. The virtual hand's postures shown on bottom side of Fig. 9 were performed using prototype V4. We see that for the last configuration, we do not have exactly the same configuration as in reality, especially for the thumb. Indeed, we have only one dof for bending and extension as a planar movement. Thus, it is not possible to reflect the adduction/abduction movement of a real thumb. To solve this problem, we have to integrate bio-mechanical based kinematics models of the hand.

For the two other cases, we can see that the virtual configurations are very close to the real ones. We added some visual markers on the virtual hand to help the user, especially when the virtual hand or a finger touches an object. Even if the vibrations sensed by the user allow him to better reach the desired configurations, adding this visual feedback increases also the user's dexterity. We measured the time needed to reach the first configuration of Fig. 9 in configurations where there is no visual and tactile feedbacks, tactile feedback only, visual feedback only and, finally, both feedbacks. Considering the virtual hand moves at 0.2m.s -1 and the object is at an initial distance of 70cm from the virtual hand, we get the following times respectively: 24s, 12s, 12s and 10s. Note that in reality, at equal displacement speed and distance, we get a time of 7s. We can clearly see the interest of having multi-sensorial feedback as, for the same task, we double the user performance between the cases without any feedback compared to both feedbacks. Indeed, without any visual or tactile feedback, it is necessary to modify the camera viewpoint to ensure a correct finger position. We observe that tactile feedback only does not bring more dexterity than just the visual one, probably because the user relies on visual markers, and we are in a configuration without any occlusion problem. However, the combination of both visual and tactile feedbacks brings advantages. Finally, Fig. 10. Evaluation of the prototype V4. The numbers along the X axis correspond to the criteria listed in Table I.

if we compare with the prototype V3s, which does not provide any tactile feedback, for the same task, we get 29s and 19s respectively. The times observed are greater than those with the prototype V4 because, as mentioned earlier, the sensors on the prototype V3s require higher forces that can generate undesired motions. Based on this comparison, we show the interest of the sensors used in the prototype V4.

Our new device significantly improves the interactions, thanks to the simultaneous use of several dofs, allowing more complex tasks with virtual avatars for example.

Finally, we also carried out the same users' tests with our new prototype. We obtained the results depicted in Fig. 10, which confirms the improvement of our prototype.

VI. CONCLUSIONS

We described and analyzed a prototype of an immersive peripheral device allowing a user to control a virtual hand and its gestures. The analysis of users' tests on previous versions of the device originated the design of this new version improving the control of the virtual hand postures as well as the dexterity to manipulate virtual objects, thanks to multi-sensorial feedbacks (tactile and visual). We especially concentrated on including tactile feedback to improve the user's immersion and interaction quality. Studies on the device shape, sensors' technologies, led to issues where physical phenomena, psycho-physiological behaviors as well as the device structure influence significantly the global interaction between real and virtual worlds. Thanks to the combination of several sensors' technologies, our prototype is able to provide relevant tactile feedback, while being comfortable. As for the previous versions of the prototype, our new version is cheap, easy to integrate in desktop environments and is calibration-free.

The next step is to set up more experiments where muscular activities of the user's arm and forearm are measured so that our prototype can be validated in a quantitative way.

The HandNavigator was designed to be integrated in several applications, such as physical simulation, industrial processes, teleoperation of robots or to communicate with virtual avatars. We will consider these applications in a near future to bring new interaction capabilities.

Fig. 1 .

 1 Fig. 1. Prototype of the HandNavigator integrating a tactile feedback.

Fig. 2 .

 2 Fig. 2. Existing prototypes (from left to right: V2, V3t, and V3s).

Fig. 3 .

 3 Fig. 3. Grasping a virtual giraffe.

Fig. 4 .

 4 Fig. 4. Evaluation of the prototypes. The numbers along the X axis are the identifiers of the criteria listed in TableI.

Fig. 5 .

 5 Fig. 5. Learning curves for each prototype in the scenario of Fig. 3. The curves do not evolve significantly after five trials.

Fig. 6 .

 6 Fig. 6. Module for one finger containing the lever-switch and pressure sensors and a vibrator.

Fig. 7 .

 7 Fig. 7. Final prototype V4.

Fig. 8 .

 8 Fig. 8. Schematic design of the final prototype V4.

Fig. 9 .

 9 Fig. 9. Different possible configurations of virtual/real hands and fingers to grasp an object.

TABLE I QUESTIONNAIRE

 I FOR EVALUATING EACH PROTOTYPE. Global motion of the virtual hand (difficult (0) -easy (

	Criteria	V2	V3t	V3s
	1.

J.-C. Léon is with Grenoble University, INRIA, Jean Kuntzmann Lab-