Mohamed Ali Kandi
email: mohamed-ali.kandi@hds.utc.fr

Hicham Lakhlef
email: hicham.lakhlef@hds.utc.fr

Abdelmadjid Bouabdallah
email: madjid.bouabdallah@hds.utc.fr

Yacine Challal
email: yacine.challal@hds.utc.fr

A Key Management Protocol for Secure Device-to-Device Communication in the Internet of Things

Keywords: Internet of Things, Device-to-Device communication, Security, Peer-to-Peer Key Management

The Internet of Things (IoT) is a network made up of a large number of devices which are able to automatically communicate in a Peer-to-Peer manner. The aim is to provide various services for the benefit of society. One of the main challenges facing the IoT is how to secure this Device-to-Device communication. Among all the security issues, the Key Management is one of the most difficult. This is mainly due to the fact that most of these devices have limited resources in terms of storage, calculation, communication and energy. Although different approaches have been proposed to deal with this problem, each of them presents its own limitations and weaknesses. In this paper, we propose a novel Key Management protocol for Device-to-Device communication in the Internet of Things. Compared to the existing Peer-to-Peer schemes, our solution provides the best compromise between the IoT requirements: resilience, connectivity, efficiency, scalability and flexibility. To achieve this balance, the network members are uniformly distributed into logical sets. A device shares then a distinct pairwise key with each member of its set and a unique pairwise set key with the members of each of the other sets. We then prove that our solution is resilient as the capture of a member compromises a negligible part of a large network. Moreover, we show that our scheme has a good network connectivity. It is then efficient as it does not require additional calculation or communication costs on the network members. We also demonstrate that our protocol is scalable as storage cost on the network members does not significantly increase when the network gets larger. We finally show that our solution is flexible.

I. INTRODUCTION

The number of devices connected to the Internet is constantly increasing since its appearance. Now that this number far exceeds that of people in the world, we are no longer talking about the Internet but about the Internet of Things (IoT). This emerging technology gives rise to revolutionary applications such as health care, environment monitoring, smart homes, smart cities, autonomous vehicles...etc. The IoT devices are able to automatically communicate to each other in a Peer-to-Peer manner. The Vehicle-to-Vehicle communication, for example, allows vehicles to exchange information about their speed and position to avoid crashes, ease traffic congestion and improve the environment. One of the main challenges facing the IoT is how to secure this Device-to-Device communication. Among all the security issues, the Key Management is one of the most difficult. This is mainly due to the fact that most of these devices have limited resources in terms of storage, calculation, communication and energy [START_REF] Alaba | Internet of Things security: A survey[END_REF].

The Key Management (KM ) is the core of secure communication. Its main role is to establish secure links between the network members [START_REF] He | Dynamic key management in wireless sensor networks: A survey[END_REF]. To achieve this, the KM provides them with secret cryptographic keys that are used to encrypt and decrypt the exchanged data. According to the encryption technique used, the KM approaches can be classified into two categories: symmetric and asymmetric. Symmetric schemes involve the use of the same key for encryption and decryption, while asymmetric approaches use two different keys. Asymmetric approaches usually imply intensive computing, which makes them impractical on the IoT constrained devices [START_REF] Yousefpoor | Dynamic key management algorithms in wireless sensor networks: A survey[END_REF].

Although symmetric schemes are more suitable for the IoT constrained devices, they rarely meets all its requirements: resilience (Rsl), connectivity (Cnt), efficiency (Efc), scalability (Scl) and flexibility (Flx). A protocol is resilient if the capture of a node does not jeopardize the communications of the other ones. When network connectivity is low, some neighboring communicators may not share a common key and relay on intermediate nodes to establish a secure link. Efficiency implies minimal use of node resources in terms of storage, calculation and communication. If increasing the network size does not significantly degrade its performance, the protocol is scalable. Finally, a flexible solution operates well regardless of nodes location and supports their dynamic deployment [START_REF] Cheikhrouhou | Secure group communication in wireless sensor networks: a survey[END_REF].

In this paper, we propose a novel KM protocol for Deviceto-Device communication in the IoT. Compared to the existing Peer-to-Peer schemes, it provides the best compromise between the IoT requirements mentioned above. To achieve this balance, the network members are uniformly distributed into logical sets. A device shares then a distinct pairwise key with each member of its set and a unique pairwise set key with the members of each of the other sets. We then prove that our solution is resilient as the capture of a member compromises a negligible part of a large network. Moreover, we show that our scheme has a good connectivity. It is then efficient as it does not require additional calculation or communication on the devices. We also prove that our protocol is scalable as storage cost does not significantly increase when the network gets larger. We finally show that our solution is flexible.

The remainder of this paper is organized as follows: related works are discussed in Section II. We detail then our solution in Section III. Section IV presents the security analysis of our protocol. In Section V, we evaluate the performance of our solution. Finally, we conclude our work in Section VI.

II. RELATED WORKS

In this work, we focus on the symmetric schemes since they require less resources than the asymmetric ones. This makes them more suitable for the IoT constrained devices [START_REF] Yousefpoor | Dynamic key management algorithms in wireless sensor networks: A survey[END_REF]. Symmetric schemes can in turn be classified into three subcategories: deterministic, pure probabilistic and deployment knowledge based schemes.

A. Deterministic schemes

The basic deterministic scheme (Pairwise Key Protocol) consists of storing a distinct pairwise key in each pair of nodes before their deployment. Since such a technique [START_REF] Choi | The best keying protocol for sensor networks[END_REF] requires a lot of storage, other approaches were proposed: Polynomialbased protocols [START_REF] Baburaj | Polynomial and multivariate mapping-based triplekey approach for secure key distribution in wireless sensor networks[END_REF][START_REF] Blundo | Perfectly-secure key distribution for dynamic conferences[END_REF], Matrix-based schemes [START_REF] Du | A pairwise key predistribution scheme for wireless sensor networks[END_REF][START_REF] Tsai | Key Management in Internet of Things via Kronecker Product[END_REF]...etc. The deterministic approaches have the advantage of being resilient against node capture. Moreover, they are usually efficient and have a good connectivity. However, they suffer from poor scalability (Figure 1a). Indeed, Pairwise Key schemes require that a node stores as many keys as there are members in the network. On the other hand, the larger is the network, the more vulnerable the Polynomial and Matrix-based approaches are to compromise. Finally, most of these schemes lack flexibility as they are based on key pre-distribution.

B. Pure probabilistic schemes

The first pure probabilistic scheme was introduced in [START_REF] Eschenauer | A key-management scheme for distributed sensor networks[END_REF]. It consists of using a large pool of keys and to randomly distribute some of them (a key ring) to each network member. Two neighboring nodes can then communicate only if they share a common key. Otherwise, they relay on intermediate nodes to establish secure links. Other methods were proposed to enhance the resilience. Using the Q-composite [START_REF] Chan | Random key predistribution schemes for sensor networks[END_REF] scheme, nodes can communicate only if they share Q keys. Also, polynomial pool based schemes [START_REF] Ruj | Pairwise and triple key distribution in wireless sensor networks with applications[END_REF][START_REF] Zhang | Key establishment scheme for wireless sensor networks based on polynomial and random key predistribution scheme[END_REF] use a pool of polynomials instead of keys. Probabilistic schemes are resilient and more scalable than the deterministic ones. Nevertheless, they suffer from poor flexibility, efficiency and connectivity (Figure 1b). Indeed, they are usually based on key predistribution. Moreover, intermediate nodes may be necessary to establish secure links. This requires additional calculation and communication and thereby more energy consumption [START_REF] Zhan | A novel key generation method for wireless sensor networks based on system of equations[END_REF]. Some works tried to enhance the connectivity using the unital design theory [START_REF] Bechkit | A highly scalable key pre-distribution scheme for wireless sensor networks[END_REF], system of equations [START_REF] Zhan | A novel key generation method for wireless sensor networks based on system of equations[END_REF]...etc. However, as long as they are probabilistic, the connectivity is rarely total.

C. Deployment knowledge based schemes

These schemes are neither deterministic nor purely probabilistic. They are based on the location of nodes to maximize the connectivity. Thus, to increase the probability of sharing keys, nodes are distributed into regional zones. Key rings are then assigned to them so that neighboring nodes share a maximum of keys. Like the other approaches, the deployment knowledge based schemes can use pairwise keys [START_REF] Choi | Location-based key management strong against insider threats in wireless sensor networks[END_REF], polynomials [START_REF] Liu | Improving key predistribution with deployment knowledge in static sensor networks[END_REF] or matrices [START_REF] Yu | A robust group-based key management scheme for wireless sensor networks[END_REF]. These approaches are resilient, scalable and provide a better network connectivity than the pure probabilistic schemes. However, they are not flexible and are more suitable for static networks (Figure 1c).

III. OUR SOLUTION

Our literature review shows that none of the existing solutions meets all the IoT requirements. To improve the scalability of deterministic schemes without loss of efficiency or connectivity, as it is the case with probabilistic schemes, our solution uniformly distributes the network members into logical sets. A node shares then a distinct pairwise key with each member of its set and a unique pairwise set key with the members of each of the other sets. The scalability of the protocol is then improved as nodes store fewer keys. Although the members of a set share the same pairwise set key, we prove that our solution remains resilient against node capture. Unlike deployment knowledge schemes, our protocol operates well regardless of the position of nodes. Moreover, as keys are dynamically distributed to the network members when nodes join or leave the network, our solution is flexible (Figure 1d).

Since some keys are shared by several nodes, the KM must ensure that they are known only by the current members. Thus, when a node joins or leaves the network, these keys are revoked and new ones are distributed to the remaining ones. This rekeying ensures that a joining node will not have access to the old keys (backward secrecy) and a leaving member will no longer know the future ones (forward secrecy).

A. Overview

The network members are distributed into logical sets to improve the protocol scalability. To each set S is associated a unique ID, sid S . It is important to note that these grouping is logical and transparent to the application layer. Although nodes belonging to the same set are considered as cognates, they can be physically far from one another. Each of them, u, is also associated with a unique ID, nid u . The keys managed by our solution can be classified into two types: Data Encryption Keys (DEKs) and Key Encryption Keys (KEKs). The DEKs are symmetric pairwise keys that are used by nodes to encrypt the data exchanged between them. The KEK are used to secure the communications between the KM and the nodes in order to protect the DEKs and thereby ensure the backward and forward secrecy. A node u, belonging to a set S , holds the following keys:

• A pairwise node key, K v u , for each of its cognates v . It is a DEK allowing the two nodes to communicate safely.

• A pairwise set key, K T S , for each set T (T = S ). This is a DEK that allows the members of S to communicate with those of T safely. This key is used, instead of the pairwise node keys, to reduce the storage overhead on nodes and thereby improve the protocol scalability.

• A node key, K u , which is a KEK known only to u. It allows the KM to communicate with the node safely. • A set key, K S , which is a KEK known to the members of S only. This key is used, instead of the node keys, when the KM sends the same message to all the members of S . The latter will then be encrypted once for more efficiency.

Figure 2 and Table I show an example of distribution of the nodes of a network N and the keys they know. Hereafter, the keys K v u and K u v are the same and can be used interchangeably. The same goes for the keys K T S and K S T . We also use the notations K R and KDF to refer to a refresh key and a pseudo-random key derivation function, respectively. Finally, KM → R :< {M } K > means that the KM sends the message {M }, encrypted using the key K , to the node(s) R.

1 K 1 K 2 1 , K 3 1 K S 1 K S 2 S 1 , K S 3 S 1 2 K 2 K 1 2 , K 3 2 3 K 3 K 1 3 , K 2 3 4 K 4 K 5 4 , K 6 4 K S 2 K S 1 S 2 , K S 3 S 2 5 K 5 K 4 5 , K 6 5 6 K 6 K 4 6 , K 5 6 7 K 7 K 8 7 , K 9 7 K S 3 K S 1 S 3 , K S 2 S 3 8 K 8 K 7 8 , K 9 8 9 K 9 K 7 9 , K 8 9

B. Rekeying upon joining

Let us consider a node u joining the network. The node is first assigned to a set S following the steps described in section III-D. The KM generates then some new keys and updates some of the previously existing ones. The aim of this update is to ensure the backward secrecy. Indeed, if these keys are not updated and if the joining node u has stored the messages previously exchanged, it will be able to decipher some of them. Next, the KM provides some nodes with the new keys and sends to others the elements allowing them to update some of the keys they hold. The process of rekeying upon joining consists of the four following steps.

1) Key generation: The first step in the rekeying process consists of determining the secret key, K u , of the joining node u. After that, the KM generates a pairwise node key, K v u , for each node v of the set S . The KM also determines the unique node ID, nid u , associated to u.

2) Key update: The KM starts by randomly generating K R . Then using it and the KDF , the KM updates the set key of S and the pairwise set keys known by its members (Formulas 1 and 2, respectively). As previously said, the aim of this update is to guarantee the backward secrecy.

K + S = KDF (K S ||K R ) (1)
K T + S = KDF (K T S ||K R ), ∀T ∈ N (2)
3) Key distribution: After the keys generation and update are completed, the KM distributes these new keys to the appropriate nodes. Thus, it sends to each node v of the set S the unicast message JM1 encrypted by means of the node secret key, K v . The message contains the ID of the joining node and the pairwise node key, K v u , associated to it. The KM also broadcasts for each set T (including S ) the message JM2 encrypted using K T , the current set key of T . The message contains the ID of the set S and K R . Finally, the KM provides u, via a pre-existing secure channel, with its secret key, the new set key, the pairwise node keys to share with its cognates and all the new pairwise set keys associated to S . After the keys distribution, the KM discards K R .

JM1 : KM → v :< {nid u , K v u } Kv > (∀v ∈ S) JM2 : KM → T :< {sid S , K R } KT > (∀T ∈ N )
4) Key installation: When a member of the set S , v , receives the messages JM1 and JM2 , it first decrypts them using its secret and set keys, respectively. Then, it installs K v u as the pairwise key to use for encrypting the communications with the joining node u. The node v also uses K R and the KDF to update the set key and all the pairwise set keys it knows (Formulas 1 and 2, respectively). After that, v discards K R . On the other hand, when a node w , not belonging to S , receives JM2 , it first decrypts the message, using the current set key, and retrieves K R . Then, using the KDF , it updates the pairwise set key it shares with the members of S (Formula 2). Once done, w discards K R .

C. Rekeying upon leaving

A node u can leave the network or be evicted when it get compromised. In both cases, the keys it knows must be revoked. The KM removes then some of them and updates some others. The aim of this update is to ensure the forward secrecy. Indeed, if these keys are not updated, the leaving node will be able to decipher some of the future communications. Next, the KM provides the network members with the elements allowing them to remove the keys that should be removed and to update those that must be updated. The process of rekeying upon leaving consists of the four following steps.

1) Key removal: The KM starts by removing the secret key, K u , of the leaving node as well as its ID, nid u . Next, the KM deletes all its pairwise keys,

K v u (v ∈ S , v = u).
2) Key update: The KM starts by randomly generating K R . Then using it and the KDF , the KM updates the set key of S and all the pairwise set keys known by its members (Formulas 1 and 2, respectively). As previously said, the aim of this update is to guarantee the forward secrecy.

3) Key distribution: After the keys removal and update are completed, the KM distributes the new keys to the appropriate nodes. Thus, it sends, to each node v of the set S , the unicast message LM1 encrypted by means of the node key, K v . The message contains the ID of the leaving node and K R . The KM also broadcasts, for each set T (T = S ), the message LM2 to provide its members with K R . The message LM2 is encrypted using K T , the current set key of T . The message LM2 is not sent to the members of S because the leaving node u knows the set key K S . The refresh key is therefore sent to the other members of S via the unicast message LM1 instead. After the keys distribution, the KM discards K R .

LM1 : KM → v :< {nid u , K R } Kv > (∀v ∈ S, v = u) LM2 : KM → T :< {sid S , K R } KT > (∀T ∈ N, T = S)
4) Key installation: When a member of S , v , receives LM1 , it first decrypts the message, using its secret key K v , and retrieves K R . Then, it removes the pairwise key K v u , which was used for encrypting the communications with the leaving member u. The node v also uses the KDF to update the set key and all the pairwise set keys it knows (Formulas 1 and 2, respectively). Once done, the node v discards K R . On the other hand, when a node w , not belonging to the set S receives LM2 , it first decrypts the message, using the current set key, and retrieves K R . Then, using it and the KDF , the node updates the pairwise set key it shares with the members of the set S (Formula 2). Finally, the node w discards K R .

D. Set management

The set management consists of distributing nodes sets while minimizing the number of keys they store. The aim is to improve the protocol scalability without significant loss of resilience, efficiency nor network connectivity. In the following, we use the notations n and p to refer to the number of nodes and sets in the network, respectively. We also denote the number of the members of a set S by m s .

A node stores one secret key, m s -1 pairwise node keys, one set key and p -1 pairwise set keys. Storage on nodes is therefore proportional to p + m s . The problem consists then of creating sets and assigning nodes to them so as to satisfy:

∀S , min (p + m s ) (3) p s=1 m s = n (4)
To have the same number of keys stored on each network member, we opted for a uniform distribution (i.e. ∀S , m s = m). By replacing 4 in 3 and studying the monotony of the resulting function (f (p) = p + n p ), we can easily show that storage is minimized when p = m = √ n. The set management aims then to uniformly distributes the n nodes of the network into √ n sets of √ n members each (Figure 2). The Assignment Algorithm (Algorithm 1) is run when nodes join the network and assigns them to the right sets. It takes as input n, the current number of network members, and assigns the joining node to a set according to the input value. The algorithm manipulates then a list of sets, ls, of size p. Each of its items contains the ID of a set and its size.

Algorithm 1: Assignment Algorithm

Input : n = the number of network members The Reorder Algorithm (Algorithm 2) is run, after a node leaving, to reduce the number of sets. It takes as input the size of the network, n, the percentage of merging, pcm, and tries to remove or merge sets when it is possible. Due to space constraints, we do not discuss the set merging in this paper.

Algorithm 2: Reorder Algorithm

Input : n = the number of network members pcm = percentage of merging

1 if m s = 0 then Remove S ; 2 else 3 if m s < pcm. √ n then 4
Find T such as m t < pcm. √ n;

5 if a set T is found then Merge S and T ; 6 end 7 end 8 Update ls;

IV. SECURITY ANALYSIS Since some keys are shared by several nodes, we start by showing that our solution fulfills the backward and forward secrecy requirements. We then prove that it provides a good level of resilience. We assume that the KM itself is secure and that only the network members can be compromised.

A. Backward and forward secrecy

The issue is to prove that a joining node cannot access the old keys and that a leaving node cannot access the new ones.

Proposition 1: Backward secrecy is guaranteed as the joining node never gets knowledge of the security material used before it joins the network.

Proof: Let us consider a node u that joins a set S of the network. The KM first updates the keys mentioned above. Then, before u can actually join the network, the KM rekeys all its current members by means of messages JM1 and JM2 . The content of the former is encrypted using the nodes secret keys and that of the latter is protected by means of the set keys. Since none of these keys are known by u, the joining node is excluded from the process of rekeying.

Proposition 2: Forward secrecy is guaranteed as the leaving node does not have access to the new security material.

Proof: Let us consider a node u that leaves a set S of the network. The KM first rekeys the cognates of u, by means of the message LM1 , then the members of the other sets, using LM2 . The content of the former encrypted using the nodes secret keys and that of the latter is protected by means of the set keys. Since none of these keys are known by u, the leaving node is excluded from the process of rekeying.

B. Resilience against node capture

According to [START_REF] He | Dynamic key management in wireless sensor networks: A survey[END_REF], resilience is the measure of the impact of one captured node on the rest of the network. The issue is then to prove that, using our solution, this impact is negligible for large networks such as the IoT.

1) Theoretical analysis: We start by analyzing the rate of compromised links due to a node capture.

Lemma 1: A node can decrypt a number of links equal to:

D = n -1 + ( √ n -1 )(n - √ n) = ( √ n -1 )(n + 1 ) (5)
Proof: A node can decrypt the communications linking it to the n -1 other network members as well as the links between its √ n -1 cognates and the n -√ n other nodes. Proposition 3: The percentage of links that a compromised node can decipher is equal to:

P = D T = ( √ n -1 )(n + 1 ) n(n-1 ) 2 = 2 (n + 1 ) ( √ n + 1 )n (6)
Proof: On the one hand, we have the lemma 1. On the other, the total number of links in a network of n nodes is equal to:

T = C 2 n = n(n -1 ) 2 (7)
Proposition 4: The compromise of the whole network occurs if and only if all the network members are captured.

Proof: Deciphering all the intra-subgroup communications requires the knowledge of all the pairwise node keys associated to it. This is only possible if all the subgroup members are captured. Also, deciphering all the inter-subgroup communications requires the knowledge of all the pairwise subgroup keys. This is only possible if at least a member of each subgroup is compromised. Deciphering all the communications is then possible if and only if all the network members are captured.

2) Comparison: According to [START_REF] Ruj | Pairwise and triple key distribution in wireless sensor networks with applications[END_REF], a perfect resilience can be reached if each pair of nodes share a distinct pairwise key. Thus, a captured node can only decipher the n -1 communications linking it to the other network members. The percentage of compromised links is equal to 2 (n-1 ) n(n-1 ) = 2 n . Due to space constraints, we only compare the resilience of our solution to the Pairwise Key schemes. Providing a perfect resilience, none of the other solutions can do better. We consider then as example the work presented in [START_REF] Choi | The best keying protocol for sensor networks[END_REF]. We proved that, using our solution, the rate of compromised links due to a node capture is equal to P . Figure 3 shows that this value is negligible for large networks such as the IoT. It is even comparable to the percentage provided by the perfectly resilient scheme. We also showed that the compromise of the whole network requires the capture of all its members. Our solution provides then a good level of resilience. Fig. 3: Variation of the rate of captured links according to n.

V. PERFORMANCE EVALUATION

Due to space constraints, we assume that the KM has plentiful of resources and focus on the protocol costs on nodes.

A. Theoretical analysis

We start by analyzing the protocol overheads on nodes. √ n keys. Moreover, regardless of the rekeying operation performed (node joining, node leaving,...etc), a node receives a constant number of messages and calculates the hash of the 2 .

√ n keys it knows.

B. Comparison

After showing that our solution provides a good level of resilience, let us prove that it meets the other IoT requirements.

1) Scalability: Although having a perfect resilience, the storage cost of the Pairwise Key schemes, in general, and the work presented in [START_REF] Choi | The best keying protocol for sensor networks[END_REF], in particular, is of the order of O(n). The deterministic solution (Kronecker) presented in [START_REF] Tsai | Key Management in Internet of Things via Kronecker Product[END_REF] and the probabilistic one (Trade) presented in [START_REF] Ruj | Pairwise and triple key distribution in wireless sensor networks with applications[END_REF] has a storage proportional to O( √ n). For the other schemes (e.g. [START_REF] Bechkit | A highly scalable key pre-distribution scheme for wireless sensor networks[END_REF]), it is difficult to deduce the storage from the network size as it depends on other parameters. Our solution also has a storage proportional to O( √ n).

Figure 4 shows that our solution stores fewer keys than the Pairwise Key schemes and can operate on larger networks of compromised nodes such as the IoT. It even provides a level of scalability comparable to the solutions presented in [START_REF] Ruj | Pairwise and triple key distribution in wireless sensor networks with applications[END_REF][START_REF] Tsai | Key Management in Internet of Things via Kronecker Product[END_REF]. Fig. 4: Variation of nodes' storage overhead according to n.

2) Connectivity: Although being scalable, the probabilistic schemes, mentioned above, suffer from poor connectivity. The probability that two neighboring nodes share a common key does not exceed 0 .25 in [START_REF] Ruj | Pairwise and triple key distribution in wireless sensor networks with applications[END_REF], while in [START_REF] Bechkit | A highly scalable key pre-distribution scheme for wireless sensor networks[END_REF] it is approximately lower bounded by 0 .632 . Using our solution, this probability is always equal to 1 . Indeed, each pair of communicators share a pairwise node or subgroup key and can establish a direct secure link. It provides then a good connectivity.

3) Efficiency: Unlike most of the probabilistic schemes [START_REF] Bechkit | A highly scalable key pre-distribution scheme for wireless sensor networks[END_REF][START_REF] Ruj | Pairwise and triple key distribution in wireless sensor networks with applications[END_REF], our solution has a good connectivity. Thus, it does not require additional calculation and communication costs to establish secure links. It also stores fewer keys on nodes than the Pairwise Key schemes [START_REF] Choi | The best keying protocol for sensor networks[END_REF] (Figure 4). Moreover, the communication and calculation costs are the same as in [START_REF] Tsai | Key Management in Internet of Things via Kronecker Product[END_REF], while they are of the order of O((log(n)) 2 ) and O(log(n)) in [START_REF] Ruj | Pairwise and triple key distribution in wireless sensor networks with applications[END_REF]. Note that communication is the operation that consumes the most nodes' energy. On top of that, our solution is based on symmetric cryptography. It is therefore efficient.

4) Flexibility: Although deployment knowledge schemes [START_REF] Choi | Location-based key management strong against insider threats in wireless sensor networks[END_REF][START_REF] Liu | Improving key predistribution with deployment knowledge in static sensor networks[END_REF][START_REF] Yu | A robust group-based key management scheme for wireless sensor networks[END_REF] provides good connectivity, they are based on nodes' location. Our solution operates well regardless of the position of nodes and supports their dynamic deployment. We previously showed that nodes can join and leave the network at any time without jeopardizing its security. It is then more flexible and suitable for dynamic networks such as the IoT.

To sum up, considering the IoT requirements (resilience, connectivity, efficiency, scalability and flexibility) as a whole, our solution provides the best compromise between them (Figure 1).

VI. CONCLUSION

In this paper, we proposed a novel Key Management protocol for Device-to-Device communication in the IoT. Compared to the existing Peer-to-Peer schemes, our solution provides the best compromise between the IoT requirements: resilience, connectivity, efficiency, scalability and flexibility. To achieve this balance, the network members are uniformly distributed into logical sets. A node shares then a distinct pairwise key with each member of its set and a unique pairwise set key with the members of each of the other sets.

We proved that our solution is resilient as the capture of a member compromises a negligible part of a large network. We then showed that it has a good connectivity. It is then efficient as it does not require additional calculation or communication. We also demonstrated that our scheme is scalable as storage on nodes does not significantly increase when the network gets larger. We finally showed that it is flexible.

In future works, we intend to decentralize the protocol in order not to have a single point of failure.

Fig. 1 :

 1 Fig.1: Symmetric Peer-to-Peer Key Management approaches.

Fig. 2 :

 2 Fig. 2: Example of a distribution of nodes.

Proposition 4 :

 4 Storage and calculation costs on nodes are of the order of O( √ n), while the communication is O(1 ). Proof: Using our solution, a node knows a secret key, √ n -1 pairwise node keys, a set key and √ n -1 pairwise set keys. It then stores in total 2 .

TABLE I

 I

: Example of keys held by nodes.

1

  Search in ls a set S such that m s <

	√	n;
	2 if no set is found then

3

Create a new set S ; 4 end 5 Assign the joining node to S ; 6 Update ls;

ACKNOWLEDGMENTS

This work was carried out and funded by INS2I STFOC project, Heudiasyc UMR CNRS 7253 and the Labex MS2T.


Bruno Lévy

Sylvain Petitjean

Nicolas Ray

Jérome Maillot

Least Squares Conformal Maps for Automatic Texture Atlas Generation

Keywords: CR Categories: I.3.3 [Computer Graphics] Picture/Image Generation, I.3.5 [Computer Graphics]: Three-Dimensional Graphics and Realism-Color, shading, shadowing and texture, I.4.3 [Image processing]: Enhancement-Geometric Correction, Texture Texture Mapping, Paint Systems, Polygonal Modeling * Alias|Wavefront

A Texture Atlas is an efficient color representation for 3D Paint Systems. The model to be textured is decomposed into charts homeomorphic to discs, each chart is parameterized, and the unfolded charts are packed in texture space. Existing texture atlas methods for triangulated surfaces suffer from several limitations, requiring them to generate a large number of small charts with simple borders. The discontinuities between the charts cause artifacts, and make it difficult to paint large areas with regular patterns.

In this paper, our main contribution is a new quasi-conformal parameterization method, based on a least-squares approximation of the Cauchy-Riemann equations. The so-defined objective function minimizes angle deformations, and we prove the following properties: the minimum is unique, independent of a similarity in texture space, independent of the resolution of the mesh and cannot generate triangle flips. The function is numerically well behaved and can therefore be very efficiently minimized. Our approach is robust, and can parameterize large charts with complex borders.

We also introduce segmentation methods to decompose the model into charts with natural shapes, and a new packing algorithm to gather them in texture space. We demonstrate our approach applied to paint both scanned and modeled data sets.

INTRODUCTION

A 3D paint system makes it possible to enhance the visual appearance of a 3D model by interactively adding details to it (colors, bump maps . . . ). If the discretization of the surface is fine enough, it is possible to directly paint its vertices [START_REF] Agrawala | 3D painting on scanned surfaces[END_REF]. However, in most cases, the desired precision for the colors is finer than the geometric details of the model. Assuming that the surface to be painted is provided with a parameterization, it is possible to use texture mapping to store colors in parameter space [START_REF] Hanrahan | Direct WYSIWYG painting and texturing on 3D shapes[END_REF]. Parametric surfaces (such as NURBS) have a natural parameterization. For other representations, such as polygonal surfaces, finding a parameterization is non-trivial. To decorate polygonal models with regular patterns, the lapped textures approach [START_REF] Praun | Lapped textures[END_REF] can be applied: local overlapping parameterizations are used to repeatedly map a small texture swatch onto a model.

A texture atlas is a more general representation (see, e.g., [START_REF] Igarashi | Adaptive unwrapping for interactive texture painting[END_REF][START_REF] Maillot | Interactive texture mapping[END_REF][START_REF] Pedersen | Decorating implicit surfaces[END_REF]). The model to be textured is partitioned into a set of parts homeomorphic to discs, referred to as charts, and each of them is provided with a parameterization. A texture atlas can be easily represented by standard file formats and displayed using standard texture mapping hardware. When used in a 3D paint system, a texture atlas should meet the following requirements:

• the chart boundaries should be chosen to minimize texture artifacts, • the sampling of texture space should be as uniform as possible,

• the atlas should make an optimal use of texture space.

The generation of a texture atlas can be decomposed into the following steps: 1. Segmentation: The model is partitioned into a set of charts. 2. Parameterization: Each chart is 'unfolded', i.e. put in correspondence with a subset of R 2 . 3. Packing: The charts are gathered in texture space.

The remainder of this section presents the existing methods for these three steps, and their limitations with respect to the requirements mentioned above. We then introduce a new texture atlas generation method, meeting these requirements by creating charts with natural shapes, thus reducing texture artifacts.

Previous Work

Segmentation into charts. In [START_REF] Krishnamurthy | Fitting smooth surfaces to dense polygon meshes[END_REF] and [START_REF] Pedersen | Decorating implicit surfaces[END_REF], the model is interactively partitioned by the user. To perform automatic segmentation, Maillot et al. [START_REF] Maillot | Interactive texture mapping[END_REF] group the facets by their normals. Several multiresolution methods [START_REF] Guskov | Normal meshes[END_REF][START_REF] Lee | MAPS: Multiresolution adaptive parameterization of surfaces[END_REF] decompose the model into charts corresponding to the simplices of the base complex. In [START_REF] Sander | Texture mapping progressive meshes[END_REF], Sander et al. use a region-growing approach to segmentation, merging charts according to both planarity and compactness criteria. All these approaches are designed to produce charts that can be treated by existing parameterization methods, which are limited to charts with convex borders. For this reason, a large number of charts is generated, which introduces many discontinuities when constructing a texture atlas.

Chart parameterization. Discrete Harmonic Map, described by Eck et al. [START_REF] Eck | Multiresolution analysis of arbitrary meshes[END_REF], are the most widely used. They are approximations of Continuous Harmonic Maps [START_REF] Eells | Another report on harmonic maps[END_REF], minimizing a metric dispersion criterion. Pinkal and Polthier [START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF] have shown the link between this criterion and another one named conformality, and have expressed both in terms of Dirichlet energy. Haker et al. [START_REF] Haker | Conformal surface parameterization for texture mapping[END_REF] describe a similar method in the specific case of a surface triangulation homeomorphic to a sphere.

The theory on graph embedding has been studied by Tutte [START_REF] Tutte | Convex representation of graphs[END_REF], where Barycentric Maps are introduced. The bijectivity of the sodefined parameterization is mathematically guaranteed. Floater [START_REF] Floater | Parametrization and smooth approximation of surface triangulations[END_REF] proposes specific weights improving the quality of the mapping, in terms of area deformations and conformality. In [START_REF] Lévy | Non-distorted texture mapping for sheared triangulated meshes[END_REF], a method is proposed to take additional constraints into account.

In all the methods mentioned above, since conformality is expressed as an indirect coupling between the parameters, boundary conditions are required, i.e. boundary nodes need to be fixed on a convex border in parameter space. Other expressions of conformality, such as the non-linear MIPS method [START_REF] Hormann | MIPS: An efficient global parametrization method[END_REF], make it possible to overcome this problem, and let the boundary nodes be free to move. However, this latter method requires a time-consuming non-linear optimization, and may get stuck in a local minima of the non-linear function. In [START_REF] Hurdal | Quasi-conformally flat mapping the human cerebellum[END_REF], Hurdal et al. propose a method based on circle packings, which are certain configurations of circles with specified pattern of tangencies known to provide a way to approximate a conformal mapping. Building circle packings is however quite expensive. The approach proposed in [START_REF] Sheffer | Param. of faceted surfaces for meshing using anglebased flattening[END_REF] consists in solving for the angles in parameter space. It results in a highly constrained optimization problem. Other methods [START_REF] Lévy | Constrained texture mapping for polygonal meshes[END_REF][START_REF] Piponi | Seamless texture mapping of subdivision surfaces by model pelting and texture blending[END_REF][START_REF] Zigelman | Texture mapping using surface flattening via multi-dimensional scaling[END_REF] can also extrapolate the border, but do not guarantee the absence of triangle flips and require interaction with the user. We introduce here a conformal mapping method, offering more guarantees, efficiency and robustness than those approaches.

In the case of texture mapping, not only the bijectivity of the parameterization should be ensured, but also its ability to make an optimum use of texture memory, and to accurately represent a signal stored in texture space. Sander et. al. [START_REF] Sander | Texture mapping progressive meshes[END_REF] describe an approach to minimize both a texture stretch criterion, and texture deviation between level of details. Since their approach is independent from the initial parameterization method, it can be applied to optimize the sampling of the parameterizations constructed by our method.

Charts packing in texture space. Finding the optimal packing of the charts in texture space is known as the bin packing problem. It has been studied by several authors, such as Milenkovic (see, e.g., [START_REF] Milenkovic | Rotational polygon containment and minimum enclosure using only robust 2D constructions[END_REF]), but the resulting algorithms take a huge amount of time since the problem is NP-complete. To speed up these computations, several heuristics have been proposed in the computer graphics community. In the case of individual triangles, such a method is described by several authors (see, e.g., [START_REF] Cigogni | A general method for recovering attributes values on simplified meshes[END_REF]). In the general case of charts, Sander et al. [START_REF] Sander | Texture mapping progressive meshes[END_REF] propose an approach to pack the minimal area bounding rectangles of the charts. In our case, since the charts can have arbitrarily shaped borders, the bounding rectangle can be far away from the boundary of the charts. Therefore, a lot of texture space can be wasted. For this reason, we propose a more accurate packing algorithm that can handle the complex charts created by our segmentation and parameterization methods.

Overview

The paper is organized as follows. Since it is our main contribution, we will start by introducing Least Squares Conformal Maps (LSCMs), a new optimization-based parameterization method with the following properties (see Section 2 and Figure 1): • Our criterion minimizes angle deformations and non-uniform scalings. It can be efficiently minimized by classical NA methods, and does not require a complex algorithm such as the ones used in [START_REF] Hurdal | Quasi-conformally flat mapping the human cerebellum[END_REF] and in [START_REF] Sheffer | Param. of faceted surfaces for meshing using anglebased flattening[END_REF]. • We prove the existence and uniqueness of the minimum of this criterion. Therefore, the solver cannot get stuck in a local minimum, in contrast with non-linear methods [START_REF] Hormann | MIPS: An efficient global parametrization method[END_REF][START_REF] Piponi | Seamless texture mapping of subdivision surfaces by model pelting and texture blending[END_REF][START_REF] Sander | Texture mapping progressive meshes[END_REF][START_REF] Zigelman | Texture mapping using surface flattening via multi-dimensional scaling[END_REF] where this property is not guaranteed. • The borders of the charts do not need to be fixed, as with most of the existing methods [START_REF] Eck | Multiresolution analysis of arbitrary meshes[END_REF][START_REF] Floater | Parametrization and smooth approximation of surface triangulations[END_REF][START_REF] Lévy | Non-distorted texture mapping for sheared triangulated meshes[END_REF]. Therefore, large charts with arbitrarily shaped borders can be parameterized. • We prove that the orientation of the triangles is preserved, which means that no triangle flip can occur. However, as in [START_REF] Sheffer | Param. of faceted surfaces for meshing using anglebased flattening[END_REF], overlaps may appear, when the boundary of the surface self-intersects in texture space. Such configurations are automatically detected, and the concerned charts are subdivided. This problem was seldom encountered in our experiments (note that as with classical methods [START_REF] Eck | Multiresolution analysis of arbitrary meshes[END_REF][START_REF] Floater | Parametrization and smooth approximation of surface triangulations[END_REF], if all the border nodes are fixed on a convex polygon, no overlap can occur). • We prove that the result is independent of the resolution of the mesh. This type of property may be usefull to reduce texture deviation when parameterizing different level of details of the same object.

In Section 3, we present a new segmentation method to decompose the model into charts. Thanks to the additional flexibility offered by LSCMs, it is possible to create large charts corresponding to meaningful geometric entities, such as biological features of characters and animals. The required number of charts is dramatically reduced, together with the artifacts caused by the discontinuities between the charts. Moreover, these large charts facilitate the use of regular patterns in a 3D paint system.

Section 4 presents our method to pack the charts in texture space. Since our segmentation method can create charts with complex borders, we pack the charts more accurately than with bounding rectangles, as in previous approaches. Our method is inspired by the strategy used by a 'Tetris' player. The paper concludes with some results, on both scanned and modeled meshes.

LEAST SQUARES CONFORMAL MAPS

In this section, we focus on the problem of parameterizing a chart homeomorphic to a disc. It will then be shown how to decompose the model into a set of charts, and how to pack these charts in texture space.

Notations

• scalars are denoted by normal characters x, y, u, v, • vectors are denoted by bold characters x = (x, y),

• complex numbers are denoted by capitals U = (u + iv),

• vectors of complex numbers are denoted by bold capitals U, • maps and matrices are denoted by cursive fonts U, X .

Figure 1: The body of the scanned horse is a test case for the robustness of the method; it is a single very large chart of 72,438 triangles, with a complex border. A: Resulting iso-parameter curves; B: The corresponding unfolded surface, where the border has been automatically extrapolated; C: These cuts make the surface equivalent to a disc; D: The parameterization is robust, and not affected by the large triangles in the circled area (caused by shadow zones appearing during the scanning process).

Conformal Maps

In this section, we quickly introduce the notion of conformal map. We will present further a new way to approximate the conformality criterion and the mathematical properties of this approximation. As shown in Figure 2, an application X mapping a (u, v) domain to a surface is said to be conformal if for each (u, v), the tangent vectors to the iso-u and iso-v curves passing through X (u, v) are orthogonal and have the same norm, which can be written as:

u v (u,v) conformal X(u,v) N iso-v iso-u
N (u, v) × ∂X ∂u (u, v) = ∂X ∂v (u, v), (1)
where N (u, v) denotes the unit normal to the surface. In other words, a conformal map is locally isotropic, i.e. maps an elementary circle of the (u, v) domain to an elementary circle of the surface.

It is possible to rewrite Equation 1 using differential operators, such as Laplace-Beltrami, as done in [START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF] and in [START_REF] Haker | Conformal surface parameterization for texture mapping[END_REF], which results in the well known cotangent weighting coefficients (see e.g. [START_REF] Eck | Multiresolution analysis of arbitrary meshes[END_REF]). The (u, v) parameters are then found to be the solution of two separate linear systems, one for u and one for v. The relation between u and v is indirectly taken into account by the right hand sides of the two systems. For this reason, this type of method requires the border to be fixed on a convex polygon. The MIPS method [START_REF] Hormann | MIPS: An efficient global parametrization method[END_REF] does not have this restriction, and expresses conformality as a relation linking the coefficients of the metric tensor. However, the resulting equations are non-linear. Another approach has been described in [START_REF] Sheffer | Param. of faceted surfaces for meshing using anglebased flattening[END_REF], based on the remark that the criterion defining a conformal mapping should be independent of a translation, rotation and scaling in parameter space (i.e. a similarity). The unknowns are the angles at the corners of the triangles. This requires a timeconsuming constrained optimization method.

Rather than discretizing the Laplace operator at the vertices of the triangulation, we instead take the dual path of considering the conformality condition on the triangles of the surface. Using the fact that a similarity can be represented by the product of complex numbers, we show how to turn the conformality problem into an unconstrained quadratic minimization problem. The u and v parameters are linked by a single global equation. This direct coupling of the u and v parameters makes it possible to efficiently parameterize large charts with complex borders, as shown in Figure 1. In this example, the cuts have been done manually (Figure 1-C), to create a large test case for the robustness of the method. (It will be shown in Section 3 how to automatically cut a model into charts homeomorphic to discs.)

Riemann's theorem states that for any surface S homeomorphic to a disc, it is possible to find a parameterization of the surface satisfying Equation 1. However, since we want to use the resulting parameterization for texture mapping, we add the constraint that the edges of the triangulation should be mapped to straight lines, and the mapping should vary linearly in each triangle. With this additional constraint, it is not always possible to satisfy the conformality condition. For this reason, we will minimize the violation of Riemann's condition in the least squares sense.

Conformality in a Triangulation

Consider now a triangulation

G = {[1 . . . n], T , (pj) 1 j n }, where [1 . . . n], n
3, corresponds to the vertices, where T is a set of n triangles represented by triples of vertices, and where pj ∈ R 3 denotes the geometric location at the vertex j. We suppose that each triangle is provided with a local orthonormal basis, where (x1, y1), (x2, y2), (x3, y3) are the coordinates of its vertices in this basis (i.e., the normal is along the z-axis). The local bases of two triangles sharing an edge are consistently oriented.

We now consider the restriction of X to a triangle T and apply the conformality criterion to the inverse map U : (x, y) → (u, v) (i.e. the coordinates of the points are given and we want their parameterization). In the local frame of the triangle, Equation 1 becomes

∂X ∂u -i ∂X ∂v = 0,
where X has been written using complex numbers, i.e. X = x+iy. By the theorem on the derivatives of inverse functions, this implies that ∂U ∂x

+ i ∂U ∂y = 0, (2)
where U = u + iv. (This is a concise formulation of the Cauchy-Riemann equations.)

Since this equation cannot in general be strictly enforced, we minimize the violation of the conformality condition in the least squares sense, which defines the criterion C:

C(T ) = T ∂U ∂x + i ∂U ∂y 2 dA = ∂U ∂x + i ∂U ∂y 2 AT ,
where AT is the area of the triangle and the notation |z| stands for the modulus of the complex number z.

Summing over the whole triangulation, the criterion to minimize is then

C(T ) = T ∈T C(T ).

Gradient in a Triangle

Our goal is now to associate with each vertex j a complex number Uj such that the Cauchy-Riemann equation is satisfied (in the least squares sense) in each triangle. To this aim, let us rewrite the criterion C(T ), assuming the mapping U varies linearly in T .

We consider a triangle {(x1, y1), (x2, y2), (x3, y3)} of R 2 , with scalars u1, u2, u3 associated with its vertices. We have:

∂u/∂x ∂u/∂y = 1 dT y2 -y3 y3 -y1 y1 -y2 x3 -x2 x1 -x3 x2 -x1 u1 u2 u3
,

where dT = (x1y2 -y1x2) + (x2y3 -y2x3) + (x3y1 -y3x1
) is twice the area of the triangle. The two components of the gradient can be gathered in a complex number:

∂u ∂x + i ∂u ∂y = i dT (W1 W2 W3) (u1 u2 u3) , where W1 = (x3 -x2) + i(y3 -y2), W2 = (x1 -x3) + i(y1 -y3), W3 = (x2 -x1) + i(y2 -y1).
The Cauchy-Riemann equation (Equation 2) can be rewritten as follows:

∂U ∂x + i ∂U ∂y = i dT (W1 W2 W3) (U1 U2 U3) = 0,
where Uj = uj + ivj.

The objective function thus reduces to

C(U = (U1, . . . , Un) ) = T ∈T C(T ), with C(T ) = 1 dT (Wj 1 ,T Wj 2 ,T Wj 3 ,T ) (Uj 1 Uj 2 Uj 3 ) 2 ,
where triangle T has vertices indexed by j1, j2, j3. (We have multiplied C(T ) by a factor of 2 to simplify the expression.)

Least Squares Conformal Maps

C(U) is quadratic in the complex numbers U1, . . . , Un, so can be written down as

C(U) = U * CU, (3)
where C is a Hermitian symmetric n × n matrix and the notation U * stands for the Hermitian (complex) conjugate of U. C is an instance of a Hermitian Gram matrix, i.e. it can be written as

C = M * M,
where M = (mij) is the sparse n × n matrix (rows are indexed by triangles, columns are indexed by vertices) whose coefficient is

mij = W j,T i √ d T i
if vertex j belongs to triangle Ti, 0 otherwise. For the optimization problem to have a non-trivial solution, some of the Ui's must be set to a priori values. Let us decompose the vector U as (U f , U p ) , where U f is the vector of free coordinates of U (the variables of the optimization problem) and Up is the vector of pinned coordinates of U, of length p (p n). Along the same lines, M can be decomposed in block matrices as

M = (Mf Mp) ,
where M f is a n × (n -p) matrix and Mp is a n × p matrix. Now, Equation 3 can be rewritten as

C(U) = U * M * MU = MU 2 = M f U f + MpUp 2 ,
where the notation v 2 stands for the inner product < v, v > (v stands for the conjugate of v).

Rewriting the objective function with only real matrices and vectors yields

C(x) = Ax -b 2 , (4) with
A = M 1 f -M 2 f M 2 f M 1 f , b = - M 1 p -M 2 p M 2 p M 1 p U 1 p U 2 p ,
where the superscripts 1 and 2 stand respectively for the real and imaginary part, v stands this time for the traditional L2-norm of a vector with real coordinates and x = (U 1 f , U 2 f ) is the vector of unknowns.

Note that A is a 2n × 2(n -p) matrix, b is a vector of R 2n and x is a vector of R 2(n-p) (the ui and vi coordinates of the vertices in parameter space that are allowed to move freely).

Properties

The above minimization problem has several fundamental properties which are proved in the appendix: • The matrix A has full rank when the number of pinned vertices, i.e. p, is larger than or equal to 2. • As a consequence, the minimization problem has a unique solution when p 2, given by x = (A A) -1 A b. The best value for p is 2, since in this case the mapping U can be fully conformal if the surface is developable (i.e. the minimum of the objective function is zero). In our experiments, we have pinned the two vertices maximizing the length of the shorted path between them (i.e. the graph diameter). • The solution to the minimization problem is invariant by a similarity in texture space. • The solution to the minimization problem is independent of the resolution of the mesh. This property is illustrated in Figure 3. • In texture space, all the triangles are consistently oriented if the pinned vertices are chosen on the boundary of T . In other words, triangle flips cannot occur.

SEGMENTATION

The segmentation algorithm decomposes the model into a set of charts. The design of the algorithm aims at meeting the following two requirements: 1. charts boundaries should be positioned in such a way that most of the discontinuities between the charts will be located in zones where they will not cause texture artifacts, 2. charts must be homeomorphic to discs, and it must be possible to parameterize them without introducing too much deformation.

For the first point, since the shading models depend on the normal, zones of high curvatures cause sharp variations of lighting. In these zones, a texture artifact will not be noticeable, since it will be negligible compared to the shading variation. Thus, to minimize artifacts, we will design the segmentation algorithm in such a way as to avoid chart boundaries in flat zones. In other words, it is suitable to generate large charts with most of their boundaries in high curvature zones.

For the second point, we will present an automatic approach, mimicking the way a user manually segments a model. Basically, the user attempts to decompose the model into parts resembling cylinders. To detect such cylinders, we will use an approach inspired by Morse theory, characterizing a function defined over the surface (see, e.g., [START_REF] Shinagawa | Surface coding based on Morse theory[END_REF]).

The next two sections present a feature detection algorithm, which finds curves corresponding to high curvature zones of the model, and a chart growing algorithm making them meet at these feature curves. Then, it will be shown how to validate the result, and subdivide the charts if needed. In the remainder of this section, we suppose that the surface is represented by a halfedge based data structure (see, e.g., [START_REF] Lienhardt | Extension of the notion of map and subdivisions of a 3D space[END_REF]).

Detect Features

The features detection phase can be outlined as follows: 1. Compute a sharpness criterion on the edges. We use here the second order differences (SOD), i.e. the angle between the normals, as in [START_REF] Hubeli | Multiresolution features extraction from unstructured meshes[END_REF]. It is also possible to use more elaborate criteria. 2. Choose a threshold τ so that a certain proportion of the edges is filtered out. In our examples, we kept 5 percent of the detected edges. 3. For each of the remaining edges, grow a feature curve by applying Algorithm 1.

Algorithm 1 attempts to anticipate the best paths, and filters out the small features caused by noise. Tagging the neighborhoods of the detected features avoids generating a large number of features in zones of high curvature. In our examples, the parameters are set

Expand Charts

Once the sharp features have been detected, the charts can be created. Our method is a greedy algorithm, expanding all the charts simultaneously from a set of seeds. It is similar to the s-source Dijkstra algorithm used in [START_REF] Eck | Multiresolution analysis of arbitrary meshes[END_REF] and to the region-growing paradigm used in computer vision. Since we want chart boundaries to meet at the level of features, the s-source algorithm is modified as follows:

• To select the set of seeds, the intuitive idea is to 'reverse engineer' the expected result. More precisely, we use the following method: a front is propagated from the borders and the feature curves detected by the previous algorithm, to compute a distance_to_f eatures function at each facet. Then, the seeds are found to be the local maxima of this distance_to_f eatures function.

• For closed surfaces without any detected feature, propagation is initialized from the two extremities of a diameter of the facets graph, as done in [START_REF] Lazarus | Level set diagrams of polyhedral objects[END_REF]. • Our s-source propagation uses -distance_to_f eatures as the priority function, rather than distance_to_seeds. The advantage of this approach is shown in Figure 4. • Charts are merged if they meet at a small distance d < ε from their seed.

In our experiments, ε = maxdist/4, where maxdist denotes the global maximum of distance_to_f eatures.

Our charts growing algorithm (Algorithm 2) uses the following data:

• distance_to_f eature is stored in each facet F , and denoted dist(F ); • for each chart C, the scalar max_dist(C) denotes the maximum distance to features for all the facets of C; • The set of edges chart_boundaries represents the borders of all charts. It makes it possible for a patch to be its own neighbor while remaining a topological disc.

As shown in Figure 5, our algorithm can detect cylindrical shapes, as the approach proposed in [START_REF] Lazarus | Level set diagrams of polyhedral objects[END_REF]. In our case, two configurations can be distinguished: • The cylinder corresponds to an extremity, such as the 'fingers' of the dinosaur's wings (Figure 5-B). The corresponding charts have the shape of a sock. This configuration is detected by com- puting the area/perimeter ratio. In this case, to facilitate the parameterization phase, a cut is added by starting from the seed and cutting the edges along the steepest-descent path. • The cylinder is non-extremal, and therefore non-capped. Algorithm 2 generates suitable boundaries, without requiring any special treatment (the resulting chart is 'rolled' around the cylinder).

Validate Charts

The so-constructed charts are then parameterized using the method presented in Section 2. After that, the following two criteria are tested:

• As mentioned in Section 2, no triangle flip can occur, and the border can be extrapolated. However, since the border can be non-convex, a new class of overlaps can be encountered. They are caused by a self-intersection of the border, as in [START_REF] Sheffer | Param. of faceted surfaces for meshing using anglebased flattening[END_REF]. Such configurations can be efficiently detected by the hardware, by drawing the parameter space in stencil mode. The stencil pixels drawn more than once correspond to overlaps. If such overlaps are detected, the corresponding chart is subdivided, by cutting it along the edges on the border of the overlapped zone, as shown in Figure 6. Note that our segmentation algorithm would not generate a single chart for the dinosaur's head (see Figure 5). In practice, the overlap problem has seldom appeared in our experiments, and was caused by tiny loops formed by the border. • Our criterion respects angles very well, as shown in the results section. As far as areas are concerned, large charts with zones of high curvature may result in large area variations. To detect these problems, the minimum and maximum model area/texture area ratio is measured over the facets. If the max/min ratio is greater than a certain threshold, the concerned chart is split, by growing two charts from the facets corresponding to the minimum and the maximum.

In the examples shown here, the threshold has been set to 2.

PACKING

Once the model is decomposed into a set of parameterized charts, it is possible to create a texture atlas by merging all the (u, v) domains of the charts. Usually, only a limited amount of texture memory is available. It is then suitable to minimize the unused space. In other words, given a set of possibly non-convex polygons, we want to find a non-overlapping placement of the polygons in such a way that the enclosing rectangle is of minimum area. The so-obtained texture coordinates are then re-scaled to fit the size of the texture.

The packing problem is known to be NP-complete (see, e.g., [START_REF] Milenkovic | Rotational polygon containment and minimum enclosure using only robust 2D constructions[END_REF] and [START_REF] Murata | Rectangle-packing-based module placement[END_REF]). Approaches based on computational geometry show good performances in terms of minimization of lost area, but are not efficient enough for large and complex data sets. For this reason, several heuristics have been proposed in computer graphics. For instance, in the method proposed by Sander et al. [START_REF] Sander | Texture mapping progressive meshes[END_REF], the bounding rectangles of the charts are packed. In our case, since the border may have an arbitrary shape, the bounding rectangle is not an accurate approximation.

For this reason, we propose a different algorithm, that packs the charts directly rather than their bounding rectangles. As in [START_REF] Azar | On 2d packing[END_REF], our algorithm is inspired by how a 'Tetris' player would operate, but without approximating the charts by their bounding boxes: 1. Each chart is rescaled to make its area in (u, v) space equal to its area in (x, y, z) space. 2. The maximum diameter of the charts are oriented vertically and sorted in decreasing order. 3. As shown in Figure 7, for each chart C, the top horizon h C (in red) and bottom horizon h ⊥ C (in pink) is computed. Rather than being represented by their bounding rectangles, the charts are approximated by the area between the two curves h C and h ⊥ C . As in classical approaches, in order to avoid unwanted blends caused by mip-mapping, an additional margin is added to the horizons. 4. The charts are inserted one by one, using the method described below.

As shown in Figure 7-A, the piecewise linear function h(u) representing the 'horizon' is maintained by the algorithm. For each chart C, the uC coordinate at the lower left corner of C is chosen in such a way that the lost space (in black) between the bottom horizon h ⊥ C of C (in pink) and the current horizon h (in blue) is minimized (see Figure 7-A). Then, the horizon h is updated using the top horizon h C of the current chart (in red). Since the parameter space is discretized into texels, it seems natural to represent the horizons by arrays of discretized values, with texture resolution. This makes the algorithm much simpler than using piecewise linear functions. For a chart C, all discrete values for uC are tested. The algorithm performs well, and takes less than one second to process all the data sets we have tested.

RESULTS

We have applied our method to different data sets, comprising meshes created with a 3D modeler (using subdivision surfaces) and scanned meshes. The results are equivalent to those obtained with MIPS [START_REF] Hormann | MIPS: An efficient global parametrization method[END_REF], but with mathematical guarantees, and can be more efficiently computed. As shown in Table 1, the stretch (see [START_REF] Sander | Texture mapping progressive meshes[END_REF]) measured on the result is of the same order as when using standard methods (e.g. [START_REF] Eck | Multiresolution analysis of arbitrary meshes[END_REF][START_REF] Floater | Parametrization and smooth approximation of surface triangulations[END_REF]). To optimize the mapping, it is easy to post-process the result of our method by the algorithm proposed in [START_REF] Sander | Texture mapping progressive meshes[END_REF]. Since the border nodes are naturally positionned (rather than arbitrarily fixed on a convex polygon), the result can be better (see Table 1).

Figure 8 shows some texture atlases. Note the presence of small charts, most of them corresponding to geometric details of the models (teeth, hoofs . . . ). This is not a problem for most paint systems, Harmonic Maps LSCM stretch (before optim.) 3.2 3.5 stretch (after optim.)

1.65 1.52 that can treat them properly. Some examples of textured models are shown in Figure 10. Table 2 shows the sizes of the data sets, the number of created charts, and the following statistics, obtained on a 1.3 GHz Pentium III (note that the timings for the packing algorithm are not included, since they are negligible):

• time to segment the model into charts;

• time to parameterize the charts. Our LSCM criterion (Equation 4) is minimized using the CG (Conjugate Gradient) algorithm. The independence to resolution suggests that a multi-grid approach would be even more efficient; • packing ratio obtained using an enclosing rectangle packing approach [START_REF] Sander | Texture mapping progressive meshes[END_REF] and our algorithm. • stretch measured before and after applying Sander et. al.'s optimization method as a post-processing (see [START_REF] Sander | Texture mapping progressive meshes[END_REF]).

The left histogram in Figure 9 shows the distribution of the angles in degrees between u and v gradient vectors. The mapping is nearly conformal in each triangle (the differences of lengthes between the u and v gradients we have measured are very near to zero). The right histogram shows the area deformations obtained with the 'Horse' data set, before stretch optimization. This histogram, showing texture area/model area ratios has been normalized, i.e. scaled in such a way that the mean value is mapped to 1. Note that since the mapping is nearly isotropic in each triangle, the L 2 and L ∞ stretch histogram (not shown here) have exactly the same appearance as the area histogram. As can be seen, even though our LSCM criterion is not designed to punish area deformations, few facets are distorted, and can easily be fixed by postprocessing using Sander et. al.'s method. The resulting texture atlases combine the advantages of LSCM (few chart discontinuities) and stretch-optimized parameterization (uniform sampling).

CONCLUSION

In this paper, we have presented a new automatic texture atlas generation method for polygonal models. Overall, we have proposed a complete and mathematically valid solution to the parameterization of complex models which proved to be more efficient and robust than existing methods and available tools in real production environments. Our segmentation algorithm, driven by detected features and inspired by Morse theory, decomposes the model into charts with natural shapes, corresponding to meaningful geometric entities. These two algorithms may have applications in other domains, such as re-meshing and data compression. We have successfully applied our technique to both scanned and synthetic data sets, making it possible to use existing 3D paint systems with them (Deep-Paint3D, Painter). In future works, to parameterize huge models, we will consider out-of-core algorithm, and analyze different numerical methods to minimize the LSCM criterion, including multigrid approaches and pre-conditioned CG. Including the stretch criterion directly into the LSCM criterion is also another possible future direction of research. We first show that the rank of M f is n -p when p 2. First note that since n n -2, min (n , n -p) = n -p if p 2 and the rank of M f is at most n -p. We assume that T is incrementally constructed with glue and join operations and prove the result by induction on the size of M f . We also assume, without loss of generality, that the p pinned vertices are concentrated in the initial triangulation. Let n i , ni -p be the dimensions of the matrix M (i) f at step i. Observe that since T is a non-degenerate triangulation, none of the coefficients Wj,T i is zero.

At step 0, the triangulation has n0 -p = 1 vertices and n 0 1 triangles. M (0) f has a single column and, since T is a proper triangulation, some of its coefficients are non-zero and it has rank 1 = n0 -p.

Assume that the property holds after step i. If step i + 1 is a join, then the number of rows of M

      M (i) f 0 . . . 0 W 1,T √ d T W 2,T √ d T 0 • • • 0 W i+1,T √ d T       .
It is now easy to see that its columns are linearly independent. Indeed, assume there are complex numbers λj such that

n i+1 j=1 λjm (i+1) j = 0, (5)
where the m (i+1) j are the column vectors of M (i+1) f

. If we look at the first n i coordinates of the column vectors, then Equation 5 reduces to n i j=1 λjm (i) j = 0, which implies that λj = 0, j = 1, . . . , ni, since M Since M f has rank n -p, both M 1 f and M 2 f have rank n -p when p 2. In turn, this implies that A has rank 2(n -p) when p 2.

A.2 Single Minimum

We now show that, when p 2, C(U) has a unique minimum.

First, notice that

∂C ∂x = 2(A Ax -A b).
Now, since the rank of the Gram matrix of A (i.e. A A) is the same as the rank of A, A A has rank 2(n -p) when p 2. Since A A is a square 2(n -p) × 2(n -p) matrix, it is thus invertible and the minimization problem has a unique solution (when p 2)

x = (A A) -1 A b.
The minimum of C(U) is zero when Ax = b, i.e. when A is invertible. Since it has full rank, this happens exactly when A is square, i.e. when n = n -p. Using the fact that n n -2, this implies that p = 2. We conclude that the mapping U is fully conformal (barring self-intersections) exactly when p = 2 and the triangulation T is built only with glue operations.

A.3 Invariance by Similarity

We now prove that if U is a solution to the minimization problem, then zU + T is also a solution, for all z ∈ C and T = (z , . . . , z ), z ∈ C. In other words, the problem is invariant by a similarity transformation.

First note that the vector H = (1, . . . , 1) is trivially in the kernel of M, since W1 + W2 + W3 = 0 in each triangle. Assume U is a solution of the problem. We get:

C(zU + T) = zz C(U) + 2zT * CU, = zz C(U) + 2z(MT) * MU = zz C(U), because T = z H is in the kernel of M. If C(U) = 0, then C(zU + T) = 0.

A.4 Independence to Resolution

We now show that if a given mesh is 'densified', then the solution to the augmented optimization problem restricted to the vertices of the initial mesh is the same. We prove this result when a single triangle T is split into three triangles, but the proof generalizes easily to a more general setting. So let v be the new vertex introduced in triangle T , i.e. as a linear combination of vertices v1, v2, v3: Assume also for the sake of simplicity that none of v, v1, v2, v3 is pinned. Call Ti (i = 1, . . . , 3) the triangle created that does not have vi as vertex. Then it is easy to see that dT i = αidT . M f is an n × (n -p) matrix. After insertion of v, the new matrix M + f is (n + 2) × (n + 1 -p). Indeed, one vertex is added, augmenting the number of columns by one, and three new triangles replace an old one, augmenting the number of rows by two. The structure of these matrices is as follows:

M f =    N f F 0 • • • 0    , M + f =       N f 0 . . . 0 L 0 P       .
where N f is an (n -1)×(n-p) matrix, F is 1×3, L is 3×3 and P is 3 × 1. If the coefficients of F are denoted by fj = Wj,T / √ dT , then it is easy to observe that the coefficients of L = (lij) and P = (pi) satisfy

lij = 1 √ αi αifj -αjfi , pi = 1 √ αi fi. (6)
The (n -p) × 1 solution to the initial problem is:

U f = (M * f M f ) (-1) M * f MpUp.
Consider the (n + 1 -p) × 1 solution to the augmented problem:

U + f = (M + f * M + f ) (-1) M + f * M + p U + p . (7)
Using the relations of Equation 6 and the fact that f1 +f2 +f3 = 0, it suffices then to observe that U + f = (U f , Uv) is the (unique) solution to 7, where Uv = α1U1 + α2U2 + α3U3.

In other words, the least squares conformal parameterization is unchanged at the old vertices and is the barycenter of the parameterizations of v1, v2 and v3 at the new vertex v.

A.5 Preserving Orientations

We now sketch the proof that least squares conformal maps preserve orientations, i.e. there are no triangle flips.

2 As a preliminary, note first that if complex numbers Wi are associated to vertices of a triangle as in Section 2.4, with vertices ordered counterclockwise, then

U 3 U V 2 V 1 V 3 V T T' V U V' U 1 u v LSCM
ζT = i W2W1 -W1W2 (8)
is positive (and equal to 2dT ). We again assume that the triangulation T is incrementally constructed with the glue and join operations. Denote the current triangulation by Ti. For the join operation, the result is trivial. We now prove the result when the current step is a glue. We use the notations of Figure 12. Let V and V be the images of T and T in parameter space. Let also Wj (resp. W j ) be complex numbers attached to T (resp. T ) and Xj (resp. X j ) be complex numbers attached to V (resp. V ). Since the local bases of two triangles of T sharing an edge are consistently oriented, both ζT -as defined in Equation 8 -and

ζ T = i W 1 W 2 -W 2 W 1
are positive. If we assume that the unfolding of Ti has no triangle flips, then we also have that ζV > 0, where ζV is defined as in Equation 8, replacing Wj by Xj. Now, writing down the equations defining U + f = (U f , Uv) as in the previous section, we find that

W 1 U1 + W 2 U2 + W v Uv = 0, (9)
where U1, U2, Uv are the parameterizations of vertices v1, v2, v (U1, U2 being unchanged by addition of v). Using the fact that W 1 + W 2 + W 3 = 0, X 1 = U2 -Uv and X 2 = Uv -U1, Equation 9 rewrites as

W 2 X 1 -W 1 X 2 = 0. ( 10
)
Using Equation 10 and the definition in Equation 8, we have:

ζ V = i X 1 X 2 -X 2 X 1 , = X 1 X 1 W 1 W 1 i W 1 W 2 -W 2 W 1 = X 1 X 1 W 1 W 1 ζ T > 0.
Thus, V is consistently oriented and the glue operation does not produce a triangle flip, proving the result.

Figure 2 :

 2 Figure 2: In a conformal map, the tangent vectors to the iso-u and to the iso-v curves are orthogonal and have the same length.

Figure 3 :

 3 Figure 3: Our LSCM parameterization is insensitive to the resolution of the mesh. The iso-parameter curves obtained on a coarse mesh (Figure A) and on a fine one (Figure B) are identical, and remain stable when the resolution varies within a mesh (circled zone in Figures C and D).

Figure 4 :

 4 Figure 4: A: Result of the features detection algorithm; B: The distance_to_seed function is not an optimal choice for driving our chart growing process; C: The distance_to_f eatures function shows iso-contours with more natural shapes; D: Result of our segmentation algorithm, driven by the distance_to_f eatures function.

  expand_feature_curve(halfedge start) vector<halfedge> detected_f eature for halfedge h ∈ { start, opposite(start) } halfedge h ← h do use depth-first search to find the string S of halfedges starting with h and such that: • two consecutive halfedges of S share a vertex • the length of S is than max_string_length • sharpness(S) ← e∈S sharpness(e) is maximum • no halfedge of S goes backward (relative to h ) • no halfedge of S is tagged as a feature neighbor h ← second item of S append h to detected_f eature while(sharpness(S) > max_string_length × τ ) end // for if (length(detected_f eature) > min_f eature_length) then tag the elements of detected_f eature as features tag the halfedges in the neighborhood of detected_f eature as feature neighbors end // if end // expand_feature_curve Algorithm 1: Features growing as follows: max_string_length = 5, which controls the size of the discontinuities to be filled, and min_f eature_length = 15.

Figure 5 :

 5 Figure 5: A: Our segmentation algorithm detects cylindrical shapes; B: An additional cut is added to 'sock-shaped' extremal cylinders.expand_charts priority_queue<halfedge> Heap sorted by dist(f acet(half edge)) set<edge> chart_boundaries initialized with all the edges of the surface // Initialize Heap foreach facet F where dist(F ) is a local maximum create a new chart with seed F add the halfedges of F to Heap end // foreach

Figure 6 :

 6 Figure 6: A: The dinosaur's head made of a single chart; B: Despite the absence of triangle flips, overlaps may occur, caused by self-intersections of the border; C: They can be removed by subdividing the chart.

Figure 7 :

 7 Figure 7: A: Our packing algorithm inserts the charts one by one, and maintains the 'horizon' (in blue) during the process. Each chart (in green) is inserted at the position minimizing the 'wasted space' (in black) between its bottom horizon (in pink) and the current horizon. The top horizon (in red) of the current chart is then used to update the horizon. B: Result on the dinosaur data set.

Figure 10 :Figure 9 :

 109 Figure 10: Hand-painted 3D models. Our LSCM method facilitates the use of procedural textures and complex patterns.

Figure 11 :

 11 Figure 11: Incremental construction of a triangulation that is topologically a disc. Left: glue two triangles along an edge. Right: join two existing vertices, creating a new triangle.

f

  grows by 1 while the number of columns is unchanged, so the rank is ni+1 -p = ni -p. If step i + 1 is a glue, a new vertex vi+1 and a new triangle T are added. Let v1 and v2 be the other vertices of T . The new matrix M (i+1) f is as follows:

  . Now the equation linking the last coordinate of the vectors m (i+1) j reduces to λn i+1 Wi+1,T / √ dT = 0, implying that λn i+1 = 0. Thus the columns of M (i+1) f are linearly independent and the matrix has full rank. The result is proved.

αi = 1 ,

 1 αi > 0.

Figure 12 :

 12 Figure 12: LSCMs preserve orientations.

Table 1 :

 1 Stretch optimization of the 'cow head' data setFigure 8: Data sets and associated texture space constructed by our method.

		dinosaur skull	bunny	horse
	vertices	14,669	16,949 34,834 48,485
	facets	14,384	15,124 69,451 96,968
	charts	43	40	23	44
	segmentation time (s)	8	17	30	43
	parameterization time (s) 10	23	95	190
	packing ratio (rectangles) 0.48	0.51	0.43	0.37
	packing ratio (our algo.)	0.55	0.55	0.6	0.58
	stretch (before optim.)	2.9	2.5	1.16	1.14
	stretch (after optim.)	1.26	1.55	1.14	1.12

Table 2 :

 2 Statistics and timings.

ACKNOWLEDGMENTS

We want to thank the Graphite development team (http://www.loria.fr/ ˜levy/graphite), especially Ben Li and Bijendra Vishal. Thanks also to the reviewers for their comments and help in improving this paper.

A PROPERTIES OF LSCMS

The minimization problem of Section 2 has several interesting properties when the number p of pinned vertices in parameter space is sufficient. In what follows, T is assumed to be homeomorphic to a disc.

A.1 Full Rank

We first show that the matrices M f and A have full rank when p 2 (p denotes the number of pinned vertices).

For this, recall that a triangulation that is topologically a disc can be incrementally constructed with only two operations (cf. Figure 11): the glue operation creates one new vertex and one new face, and the join operation creates one new face. Thus, incremental construction creates at most as much vertices as faces. Since the simplest triangulation (one triangle) has one face and three vertices, we have that n n -2 (where n denotes the number of vertices, and n the number of triangles, as in the rest of the paper).