Trinh Dang

Nguyen

Bao Quoc

Duong

Eric Zamai

Muhammad Kashif
email: kashif.shahzad@grenoble-inp.fr

Trinh Nguyen

Quoc-Bao Duong

Muhammad

Dynamic structure identification of Bayesian network model for fault diagnosis of FMS

Keywords: Fault diagnosis, Flexible Manufacturing Systems, Logical diagnosis, Bayesian network

This paper proposes an approach to accurately localize the origin of product quality drifts, in a flexible manufacturing system (FMS). The logical diagnosis model is used to reduce the search space of suspected equipment in the production flow; however, it does not help in accurately localizing the faulty equipment. In the proposed approach, we model this reduced search space as a Bayesian network that uses historical data to compute conditional probabilities for each suspected equipment. This approach helps in making accurate decisions on localizing the cause for product quality drifts as either one of the equipment in production flow or product itself.

I. INTRODUCTION

The manufacturing processes evolve to adapt with the increasing demand diversities. This has resulted in complex manufacturing environment where it is difficult to control, diagnose and recover the abnormal events in timely manner. The fault detection and diagnosis techniques have been a special focus in the manufacturing industries to improve product quality and reduce associated costs. Therefore, different methods are developed for more accurate and quick fault diagnosis to improve product quality such as [START_REF] Chanthery | Monitoring and active diagnosis for discrete-event systems[END_REF], [START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF]. An advantage of these methods is the flexibility to adapt to different production systems.

However, modern production systems comprise a lot of equipment that are monitored by thousands of sensors. This result in huge amount of information and makes it difficult to identify the model for fault diagnosis with corresponding variables. The artificial intelligence methods are practical techniques capable of using large amount of information from large production databases as presented in [START_REF] Verron | Fault diagnosis with bayesian networks: Application to the tennessee eastman process[END_REF], [START_REF] Weber | Overview on Bayesian networks applications for dependability, risk analysis and maintenance areas[END_REF]. These approaches can be performed without understanding the underlying structure of a production system [START_REF] Bouaziz | Dependability of complex semiconductor systems: Learning Bayesian networks for decision support[END_REF]. Especially, the Bayesian network (BN) approaches are widely used to compute the conditional probabilities between a detected fault and a set of potential failure candidates (root causes). This set and these relations between variables are normally modeled through a graphical structure during learning phase of BN. The associated probabilities help to make accurate and informed decisions about corrective maintenance. However, due to the complexity of present-day manufacturing systems, this graphical structure identification is a complex task to be performed by a personal operator [START_REF] Bouaziz | Dependability of complex semiconductor systems: Learning Bayesian networks for decision support[END_REF]. In addition, it also depends on the database size and the experts knowledge [START_REF] Freitas | Data mining and knowledge discovery with evolutionary algorithms[END_REF]. Besides, the diversity of system characteristics (structure, product type, recipe and the human factor) and the dynamic feature of production systems add to the difficulty in BN structure identification. In fact, the system elements are frequently changed due to the introduction of new manufacturing technologies and management philosophies. Due to a change in one element may trigger effect to others, we must retrain the learning step of a BN model to adapt with newly emerging situations in a production system. In this case, the time and work-load for control system are very large.

In summary, the paper focuses on an approach to dynamically generate the structure of BN and associated probabilities from the reduced search space of suspected production equipment as identified by the logical diagnosis model. For this context, we propose the use of logical diagnosis model [START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF] to limit the space of root causes and provide the cause-consequence relationships for the graphical structure identification. Thus, it is followed by learning conditional probabilities in resulting BN to accurately identify the cause of product quality drift as either a suspected equipment in the production flow or product itself. This approach enables real time diagnosis and subsequent corrective maintenance in a flexible manufacturing system (FMS).

The next section will introduce the characteristic of the case study and objectives for real time diagnosis in FMS. Section 3 presents state of the art fault diagnosis models being used in complex manufacturing environments. The proposed diagnosis presented and discussed in Section 4. The proof of concept is presented in Section 5. We conclude this paper with conclusions and future perspectives in Section 6.

II. FLEXIBLE MANUFACTURING SYSTEM

In the manufacturing industry, speed, reliability, flexibility, cost, rapid product innovation and quality are all related to design a manufacturing process [START_REF] Brown | Strategic manufacturing for competitive advantage: transforming operations from shop floor to strategy[END_REF]. Many types of manufacturing system as presented in [START_REF] Phillips | Manufacturing Plant Layout: Fundamentals and fine points of optimum facility design[END_REF] support the design and management of manufacturing processes. Following the variety of the customer demands, new products frequently are introduced. This requires FMS [START_REF] Papadopoulos | Analysis and design of discrete part production lines[END_REF] that is widely used in large scale production industries like semiconductor and automotive industries. A FMS is a complex automated production system which consists of several production stations each specializing in particular operations with an appropriate transport system. These production stations and transport system (controlled system) are controlled by a control system. The computer system coordinates the activities and executes the demands based on control laws. The flexibility of FMS is generally demonstrated by multiple products, multiple lines, multiple recipes, and multiple human factors. Generally, a complex system may have many production processes. The production process is loaded corresponding to each product and executed by the control system. The structure of a production process describes type of product, product lines and corresponding equipment in operating part of controlled system. Any change in the existing product or product line, requires customization of existing process.

Our research is focused on the general framework of a complex system with dynamic production environments. This case study has the characteristics: multiple products, multiple lines, multiple recipes, and multiple human factors as illustrated in Fig. 1. In complex engineering applications, systems are composed of many components and subsystems, and the ways in which these elements interact will affect the way failures propagate within subsystems and across subsystem boundaries [START_REF] Hine | A review of diagnostic techniques for ISHM applications[END_REF]. For monitoring the execution of these elements, the hierarchical and modular control are often used as presented in [START_REF] Jones | A multi-leve/multi-layer architecture for intelligent shopfloor control[END_REF]. In this context, an automated manufacturing system is organized by Computer Integrated Manufacturing (CIM) architecture that contains : controlled system, product flow and control system as shown in Fig. 2. The controlled system con-sists of actuators and sensors. These sensors allow controlling and monitoring the executions of actuators and product flows. Each elementary part of the controlled system (actuators and sensors) is controlled by a local control module, and these two elements are called Functional Chain (FC) [START_REF] Henry | Controlled system model adapted for control law synthesis[END_REF]. These FCs receive and execute demands from the coordination level of control system. When a FC cannot correctly execute a demand, it implies that a fault is produced. A fault may be created by a component and can propagate to other components through product lines. Therefore, this failure propagation may make consequences on many components; and so, the new faults may continue to occur. Thus, when a fault is detected by a metrology, its root causes may come from one or a part of elements in the production system. In short, it is critical to precisely and quickly locate the causes of detected fault to save recovery time (for return to a normal status) of the production system.

III. GENERAL SITUATION

For fault diagnosis in a production equipment or a simple system, the diagnosis methods such as Fault Detection and Isolation (FDI) [START_REF] White | Automated diagnosis of feature model configurations[END_REF], and Failure Mode and Effects Analysis (FMEA) [START_REF] Rao | Dynamic fault tree analysis using Monte Carlo simulation in probabilistic safety assessment[END_REF] are classically used. An advantage of these methods is their simple integration in the production system. These approaches generally model all possible failures in a production system using detection information and failure symptoms, collected through sensors. However, due to the complexity of a manufacturing system, it is too difficult to describe all the failures and so their causes. When an unpredicted fault occurs, the diagnosis cannot provide a result without expert support. In fact, a fault can propagate from an equipment to others equipment through product flows and it may have the consequences to others system elements. The unpredicted failures occur quite often; hence, the existing failure model cannot explain such events.

Following failure propagation analysis, diagnosis approaches have been proposed such as propagation graph [START_REF] Abdelwahed | System Diagnosis using Hybrid Failure Propagation Graphs[END_REF] and temporal chronic [START_REF] Strasser | Diagnostic alarm sequence maturation in timed failure propagation graphs[END_REF]. These approaches are based on historical production data to locate operations at the possible origins of detected fault. However, the problems of these approaches are demonstrated as follow: First, the actions of system operations in [START_REF] Abdelwahed | System Diagnosis using Hybrid Failure Propagation Graphs[END_REF] are considered as independent. In fact, several operations may use the same functional chain to execute the demands of control system. Indeed, the corresponding actions are not independent. Second, these approaches do not analyze the behaviors of the FCs, so they cannot explain the consequences between elements of a FC and between different FCs. And finally, these approaches cannot reduce the size of the model. Due to the cyclical operation of the control system, a large amount of information of this production system provokes the problem of combinatory explosion.

Consequently, we are especially interested in the logical diagnosis model proposed in [START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF]. In this model, a diagnosis function is proposed to characterize historical production data of a controlled system to search suspected origins against the detected fault, in real time. Thus, this model provides a set of possible origins, and its size is reduced by analyzing sensors signals from the controlled system. This is appropriate to diagnose the faults of production equipment and products. However, this model does not show suspect level of each member in this set of fault origins. It is difficult in deciding a maintenance order. In a complex system, the set of possible fault origins is still large and this method would be extended to optimize corrective maintenance activities.

To evaluate the suspect level of fault origins, the probabilistic approaches are widely used such as Neural network [START_REF] Khomfoi | Fault Diagnostic System for a Multilevel Inverter Using a Neural Network[END_REF] and Bayesian network [START_REF] Verron | Fault diagnosis with bayesian networks: Application to the tennessee eastman process[END_REF]. These probabilistic approaches have capacity to provide the probability values from the large database and the number of corresponding variables of a production system. These probability values allow evaluating the suspect level (high or low) in order to support decision of a maintenance strategy. In particular, the BN models have the advantages which fit to be applied in manufacturing industry as explained in [START_REF] Weber | Overview on Bayesian networks applications for dependability, risk analysis and maintenance areas[END_REF]. The methods based on BN are introduced in [START_REF] Bouaziz | Dependability of complex semiconductor systems: Learning Bayesian networks for decision support[END_REF] to make a diagnosis for a multiple variables system. The Confidence level of feedback information is proposed in [START_REF] Duong | Confidence estimation of feedback information for logicdiagnosis[END_REF] to provide the probability value which shows the correct actions confidence of reported information of equipment. When the database is available, these approaches are strong tools for the corrective maintenance.

However, the probabilistic approaches must be extended. First, the structure of model for Confidence level of feedback information in [START_REF] Duong | Confidence estimation of feedback information for logicdiagnosis[END_REF], is static with seven parameters. In the context of flexible manufacturing system with the characteristics such as multiple products, multiple lines, multiple recipes and multiple human factors, the production environment often change as we presented in above section. Thus, the set of parameters which can have impact on the equipment confidence is dynamic. And second, the BN model must be rebuilt with new contextual information from the production database. In practice, when the database is available, learning approaches are often used for modelling BN as presented in [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. In the learning approaches, a graphical structure and probabilistic rules are estimated from observed data. We can cite a number of learning approaches introduced in [START_REF] Freitas | Data mining and knowledge discovery with evolutionary algorithms[END_REF], [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] such as EM, Maximum likelihood, IC&IC, K2, and Genetic algorithm. Many studies in [START_REF] Bouaziz | Dependability of complex semiconductor systems: Learning Bayesian networks for decision support[END_REF], [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] show that these learning approaches are still complex in identification the variables for a human agent because they depend strongly on expert opinions. Furthermore, the learning workload for a calculating operator is still large [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. This spends too much time and is not appropriate to on-line structural identification due to this depends on the exploitation of databases, the experiences of experts and the capacity of a calculating operator. In fact, production environments are increasingly stressed by strong competition. It shows that the time for locating the root causes of failures and process recovery (return the process to a normal status) is very important. These challenges promote the researches to apply BN model for real time fault diagnosis and corrective maintenance optimization. The proposed model consists of: (i) searching possible root causes against a detected fault in the past evolution of the operating part in controlled system, (ii) computing the probability values that show the suspect level against possible root cause. This diagnosis model operates as follows: First, this model is generated by the coordination level of a control system. This coordination level sends commands and receives real time reports for and from all system components. This coordination level also provides the necessary information for the logical diagnosis. Second, once a failure is detected by a metrology, a set of possible origins and its correlations are defined by the logical diagnosis model. Third, the information of this set is sent to BN model. Thus, a graphical structure of failure mode is determined to support the structural identification in learning step of BN model. Finally, based on the theories of BN, the corresponding conditional probabilities of all nodes in network are calculated from historical production data. All calculated results are memorized in the database.

The model for diagnosis is based on the control system and dynamic reconfiguration as proposed in [START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF], [START_REF] Henry | Controlled system model adapted for control law synthesis[END_REF]. The reactive loop is characterized by collaboration of several supervision, monitoring and control (SM&C) functions such as detection, diagnosis, prognosis, decision and automatic control [START_REF] Zamai | An architecture for control and monitoring of discrete events systems[END_REF]. The coordination level of control system has capacity to manage a set of Functional Chains and corresponding information. It also can demand and memorize the information on the database of system. When a fault is detected, this model will demand the diagnosis results from logical diagnosis model and the BN model. We will explain its executions below.

The Logical Diagnosis model proposed in [START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF] is one part in the treatment of failure propagation through a product in a complex production process. In this model, a diagnosis function is proposed to characterize the historical information of a controlled system in locating the possible origins of a detected fault. This diagnosis function is presented in three main points: First, a model generated real time in normal operation and collects all the necessary information from coordination level for fault diagnosis through the Operation Models which are presented in [START_REF] Henry | Controlled system model adapted for control law synthesis[END_REF]. An Operation Model consists of information of Functional Chains (equipment, sensors and local control modular), (Pre-)conditions, (Pre-)constraints, effects of operations. This model describes a graphical structure of a production process which consists of system components and its relationships following product flows. Second, a mechanism to reduce this model is developed based on the exploitation of controlled system observations. Following an operation in this model, if the information provided by these two elements is not coherent, the controlled system sends a faulty execution report to the coordination level. If the coordination level does not receive a faulty execution report from these Functional Chains, it can conclude that the corresponding element is reliable. These reliable elements will be removed from the model, while the suspect elements will be retained. And finally, a mechanism is defined, based on the failure propagation approach, that allows us searching for the possible origins and the possible other consequences of a fault when the model receives a faulty execution report. Consequently, this logical diagnosis model provides a set and a reduced structural model of suspect operations which have logical relationships with a faulty execution. These suspect operations are considered as possible origins {O 1 , ...O i , ...O n } of detected fault. The reduced model describes the logical links between possible fault origins. The doubt can be propagated from a possible origin to another one. Hence, these logical links are considered as the cause-consequence relationships. This set of possible fault origins and the information of reduced model will be sent to BN model, and so, they are used for structural identification. Therefore the conditional probabilities are computed next.

The BN model is based on the BN theories as presented in [START_REF] Weber | Overview on Bayesian networks applications for dependability, risk analysis and maintenance areas[END_REF]. When this model receives a set of possible fault origins and reduced structural model through the coordination level, it will be performed as follows:

• First, an algorithm as presented below is used for equivalent transition from a set of possible root causes • Third, after a graphical structure is identified, the learning step of BN model is performed to calculate the probabilities P (O i | P roduct) of each node O i following the product, and next the conditional probabilities between the child-nodes and its parents based on the historical information of production system from the database. We consider that each member of the set of fault origins has two states {1, 0}. So, the conditional probabilities over each operation O j with (j = 1 as shown in equation ( 1):

for i = 1 to n do Initialize k i = 1 for i ′ = 1 to n and i ′ = i do if O i ′ is located in front of O i and has directly logical relationship with O i then O i ′ =Parent of O j O ki i = O i ′ and k i = k i + 1 end if end for N i = k i and delete k i end for return P arents(O i )= O ki i ; k i = 1 to N i ; i = 1 to n • Second,
               p O j = 1 | O 1 j = 1, ...O kj j = 1, ...O Nj j = 1 . . . p O j = 1 | O 1 j = 0, ...O kj j = 0, ...O Nj j = 0 p O j = 0 | O 1 j = 1, ...O kj j = 1, ...O Nj j = 1 . . . p O j = 0 | O 1 j = 0, ...O kj j = 0, ...O Nj j = 0                (1)
• In this paper, in order to compute the conditional probabilities, we need to extract the distribution over some subset of variables or a single variable, and thus, we need to marginalize or sum out the variables other than the variables of interest as explained in [START_REF] Manfredotti | Modeling and inference with relational dynamic Bayesian networks[END_REF]. The marginalization rule for any sets of variable X and Y is given by :

P (X) = y P (X, y); y ∈ Y (2)
The distribution over X can be obtained by summing out all the other variables from any joint distribution containing X. We can use the conditional probabilities instead of joint probabilities to compute the probabilities over X as shown in equation (3): Hence, the previous sections show that a set of possible origins is dynamically determined when a fault is detected. And when the database is available, the corresponding conditional probabilities are calculated by diagnosis model. This set of possible origins and these probability values are sent and memorized on the database by the coordination level. They are used to support decision in maintenance operations.

P (X) = y P (
The advantages of proposed model are: first to locate in real time a possible fault origins set, second to reduce the space of the possible fault origins set by the evaluation of suspect level, and finally less workload for structure identification of BN model. In this model, the set of possible root causes is significantly reduced by a logical diagnosis model. Next, this is used to simplify a structural identification of BN model. Thus, the BN model receives a graphical structure with only elements at the possible origins of a detected fault. Other elements not related to the detected fault are rejected. Consequently, it does not need to compute all probabilities corresponding with all elements of system. Therefore, the probabilistic values given by BN model help us to continue reducing the set of possible origins. This implies that we can save the recovery time of a production system. This also shows that the combination of a deterministic approach (logical diagnosis) and a probabilistic approach (BN) help us to locate more precisely and quickly equipment at that the root causes of a detected fault. Consequently, this proposed model is appropriate to apply for fault diagnosis in a flexible manufacturing system with large space of the database and identification a solution to the more approaching the optimal model. In this paper, we propose an example to explain the application of proposed diagnosis model in the case of a production system in Section 5.

V. APPLICATION EXAMPLE

In this section, an example of a production system as described in Fig. 5 consists of 4 machines and one metrology. We consider that this production system is integrated with a diagnosis model which has a diagnosis mechanism as proposed in above section. A sensors system is installed to detect the quality of each machine (good or bad) and send it to the coordination level of a control system. The product quality (good or bad) is detected by metrology. If the quality of a product is bad, it concludes that this product has a fault. In this example, we consider that the failure can propagate from one machine to another one following the product flows. Hence, we will locate the machines at the possible root causes of this detected fault. The logical diagnosis model determines an evolution chain through the Operation Models of all machines. Hence, a model of machines evolutions is constructed as shown in Fig. 6. The reliable elements (black nodes) are removed, while the suspect elements (white nodes) are retained based on the exploitation of sensors system. These describe a directed graphical model corresponding with the product flow. When a fault is detected at the metrology, following the back failure propagation, we obtain a set {M 1, M 2, M 4} that are the possible fault origins. Where:

• P(M1=bad|Product), P(M2=bad|Product) are the probabilities of M1 and M2 qualified as bad corresponding to the product.

• P(M4|M1, M2) is the conditional probability of M4 in bad state given by its parents M1, M2.

• P(fault|M4) is the conditional probability of product in bad state given by its parent M4.

Consequently, the conditional probabilities of machines (given when product is in bad state) are computed as shown in TABLE I. The result in Fig. 7 shows that the possible origins of the fault detected on product by metrology may be one machine of the set {M1, M2, M4}. Furthermore, the result shown in TABLE I concludes that the conditional probability P(M4=bad|fault) of node M4 is maximum value. Consequently, the machine M4 must be checked in the first place by a maintenance operator. This example also implies that if we do not use the BN, the set of possible origins have 3 members {M1, M2, M4}. Thus, we must detect all these 3 members. While if the BN model is used, the computed conditional probabilities helps to make decision what the member has highest suspect level. In addition, if the diagnosis process is performed based on a BN model without the Logical Diagnosis model, all the machines in Fig. 5 may be possible origins of detected fault. We must use a learning method to identify the graphical structure of the BN. The computation time and work-load become larger than the case of within Logical Diagnosis model. Indeed, the combination of two techniques allows making use of its advantages and avoiding inconveniences. The proposed model is dedicated to advance in real time diagnosis and simplify the decision making for a maintenance strategy. It is feasible to apply on the flexible manufacturing systems under dynamical production environment.

VI. CONCLUSION

A diagnosis model for accurate fault diagnosis and corrective maintenance in a complex manufacturing system has been proposed. This model is based on a logical diagnosis approach which can limit a part of equipment at the possible origins of detected fault. The logical diagnosis approach is extended with a BN model to evaluate the suspect level of each member of a set of possible fault origins. Besides, the logical diagnosis also helps to simplify the failure model identification in the learning phase of BN.

The proposed diagnosis model is appropriate in a complex manufacturing environment which has the characteristics such as: multiple products, multiple lines, multiple recipes, multiple human factors. The elements of system often change under flexible production activities. The proposed model has the advantage of dynamic localization of potential root causes, in less time and less workload is needed to compute the conditional probability values.

This study is first based on precisely and quickly locating the faults which can come from equipment and products. Moreover, the approach must take into account the different fault sources as recipes and human factors. Consequently, our future work is focused on the fault diagnosis of a general manufacturing system with multiple-fault sources.

Fig. 1 .

 1 Fig. 1. A Flexible Manufacturing System

Fig. 2 .

 2 Fig. 2. CIM architecture of an automated manufacturing system

  In this paper, we propose a diagnosis model that allows real time localization of the possible fault root causes, thus dynamically calculating the conditional probabilities between a fault and its causes. This model is based on a logical diagnosis model and a BN model. It has four main points: First, the logical diagnosis model provides a set of possible fault origins. The relationships of members in this set are used to construct a graphical structure. Second, we use this graphical structure for the structural modeling of BN model. This idea is to simplify the variables identification in learning step of the BN. Third, the historic information from the production system is used to estimate probabilistic rules in learning step of BN model. And finally, the conditional probabilities of nodes in structural model are calculated. These probabilistic values allow evaluating suspect level of each possible fault origin. We will explain in detail this diagnosis model in the next section. IV. DIAGNOSIS MODEL This paper presents a diagnosis model which uses the results given by a logical diagnosis model to dynamically identify the BN structure for fault diagnosis in a flexible manufacturing system. This model comprises a logical diagnosis model and a BN model as shown in Fig.3.

Fig. 3 .

 3 Fig. 3. Diagnosis model

  we obtain a graphical structure of BN model from previous step. Each node in this structure may be a parent of child-nodes and may be a child of other parent-nodes. For instance, a set of parents a detected fault node is a finite number of mutually exclusive operation {O 1 , ...O j , ...O N }. Each operation O j with (j = 1 • • • N ) has a set of parents as O 1 j , ...O kj j , ...O Nj j , with k j = 1 • • • N j . Therefore, each member O kj j has a set of parents as O j with l j = 1 • • • L j . In this case, the BN model has a hierarchical structure as shown in Fig.4.

Fig. 4 .

 4 Fig. 4. Graphical structure of BN model

  X|y).P (y); y ∈ Y (3) Finally, the model computes the conditional probabilities P (O i | f ault, P roduct) over node O i with (i = 1 • • • n) given by detected fault according to the Product as illustrated in the next equation: P (O i | f ault, P roduct) = n j=1 Oj P (O i ).P (O j | O i ). P (f ault | O 1 , ..., O j , ..., O n , P roduct) (4) with j = 1 • • • n, j = i.

Fig. 5 .

 5 Fig. 5. A production system

Fig. 6 .

 6 Fig. 6. Machines evolutions model following the product flows

Fig. 7 .

 7 Fig. 7. The graphical structure and the conditional probabilities of its nodes

  THE PROBABILITY RESULTS P (M 1 = bad | f ault) P (M 2 = bad | f ault) P (M 4 = bad | f ault)

  and its graphical structure given by the logical diagnosis model to a graphical structure of BN model. In Algorithm 1, we assume that a set of possible root causes has n members as {O 1 , ...O i , ...O n } obtained by logical diagnosis model. Each member O i of the set is considered a node in the BN. The other membersO i ′ ; (i ′ = 1 • • • n, i ′ = i)which are in front of the node O i following a product flow and have directly logical relationship with it, are parents of O i .

	Algorithm 1 Equivalent transition of a graphical structure
	Inputs : A set of n possible origins and its relationships.
	Outputs : A graphical structure of BN model.

  • • • N ) and its parents O 1 j , ...O

	kj j , ...O	Nj j
	are defined as a matrix P O j | O 1 j , ...O	kj j , ...O	Nj j

TABLE I .

 I


Martin Guay

Marie-Paule Cani

Rémi Ronfard

The Line of Action: an Intuitive Interface for Expressive Character Posing

Keywords: CR Categories: I.3.6 [Methodology and Techniques]: Interaction techniques-[I.3.7]: Computer Graphics-Animation Line of action, sketch-based animation, character animation

Figure 1: Expressive character poses created in a few seconds each, by sketching intuitive lines of action.

Introduction

Because humans have been drawing and sketching for centuries, many researchers-from within and from outside the computer * martin.guay@inria.fr graphics community-have proposed and argued for sketching as an intuitive and natural interface for both modeling and editing 3D virtual objects. By designing these interfaces closer to the creative process involved when sketching, the cognitive workload, i.e. the mental steps required to achieve a specific task, can be efficiently reduced. Looking at these stages-i.e. the ones involved in the creative sketching process [START_REF] Goldschmidt | The dialectics of sketching[END_REF][START_REF] Cherlin | Sketch-based modeling with few strokes[END_REF]]we usually find laid down in early stages, principal descriptive primitives-circles, lines, ovals, etc.-that give the overall perspective, shape and mass of objects and characters. For modeling 3D objects, primitive shapes such as cylinders [START_REF] Gingold | Structured annotations for 2D-to-3D modeling[END_REF] and abstract forms such as contours [START_REF] Igarashi | Teddy: a sketching interface for 3d freeform design[END_REF][START_REF] Karpenko | Smoothsketch: 3d free-form shapes from complex sketches[END_REF] have been successfully used to ease the sketching process by making coarse shape design easier.

When it comes to character posing however-a task commonly done by manipulating an articulated skeleton used to parametrize the character's geometry-sketch-based modeling research has mostly relied on stick figures as the 2D representation of the character's skeletal parametrization [START_REF] Davis | A sketching interface for articulated figure animation[END_REF]]. Although drawing 2D stick figures accelerates posing, it requires sketching multiple, intersecting strokes for specifying the individual limbs, and getting expressive poses requires skills.

Inspired by the practice of cartoonists [START_REF] Lee | How To Draw Comics The Marvel Way[END_REF][START_REF] Blair | Cartoon Animation[END_REF][START_REF] Hart | How to draw animation[END_REF]; [START_REF] Brooks | The Complete Cartooning Course[END_REF], we explore the alternative of a single, smooth stroke-the line of action (LOA)as an abstraction of the character's body. We claim that the LOA is an intuitive interface that helps creating more expressive poses. Controlling the whole body with a single stroke makes early design easier and less time consuming. And the resulting poses are more expressive as they often exhibit an aesthetic curved shape that conveys the full body expression more clearly in a given viewpoint. Examples of expressive poses rapidly created with lines of action are shown in Fig. 1. They mimic the poses of a fashion model, of two dancers, a full body swing and a super hero on the run.

Since the line of action is an abstraction of the body and the correspondence between both is not explicitly given, posing a 3D virtual character from an arbitrary LOA is a challenge. In this paper, we give a formal definition of the LOA, enabling us to solve character posing by solving an optimization problem. Our method includes an automatic way of determining the correspondence between the line of action and a subset of the character's bones. We then propose extensions to the LOA concept for drawing secondary lines and we address the well known depth ambiguities in a new, but simple way; by constraining the transformations to the viewing plane. We validate this approach by generating both expressive static poses and quickly drafting keyframe animations.

Related Work

Standard character posing methods: Posing virtual characters by manipulating a skeletal parametrization of the character's geometry is common practice in computer animation with previous works dating back to [START_REF] Burtnyk | Interactive skeleton techniques for enhancing motion dynamics in key frame animation[END_REF]. Positions on the kinematic chains are non linear and manipulating every joint is cumbersome. Inverse kinematics (IK) [START_REF] Zhao | Inverse kinematics positioning using nonlinear programming for highly articulated figures[END_REF][START_REF] Girard | Computational modeling for the computer animation of legged figures[END_REF], allows the user to focus on end-effector position targets, while the system automatically solves for the degrees of freedom in the kinematic chains. However, manipulating 3D widgets can be confusing and time consuming. For this reason, several recent works propose to sketch the characters in 2D.

Sketch-based character posing: Unlike sketch-based freeform shape modeling-which holds several high-level shape descriptors such as visible contours [START_REF] Karpenko | Smoothsketch: 3d free-form shapes from complex sketches[END_REF][START_REF] Igarashi | Teddy: a sketching interface for 3d freeform design[END_REF] or annotated shape primitives [START_REF] Gingold | Structured annotations for 2D-to-3D modeling[END_REF]]-sketchbased character posing has relied mostly on stick figures [START_REF] Davis | A sketching interface for articulated figure animation[END_REF][START_REF] Mao | A sketch-based gesture interface for rough 3d stick figure animation[END_REF][START_REF] Lin | A sketching interface for sitting-pose design[END_REF][START_REF] Wei | Intuitive interactive human character posing with millions of example poses[END_REF][START_REF] Choi | Retrieval and visualization of human motion data via stick figures[END_REF]. In contrast, our work makes use of a higher-level abstraction-the line of action-enabling us to shape characters using one or two smooth strokes, each drawn in a single hand gesture. This makes the posing interface closer to the early stages of the creative process where coarse strokes are used to describe the principal body shapes [START_REF] Goldschmidt | The dialectics of sketching[END_REF]]. It also increases consistency as the lower and upper bodies typically fit onto the same aesthetic curved stroke.

Re-constructing a 3D pose from a 2D sketch is, due to depth, highly ambiguous and under-constrained; many poses exist for the same sketch. A first solution is to let the user choose from several possible solutions as [START_REF] Davis | A sketching interface for articulated figure animation[END_REF] suggests. Another approach is to supplement the problem with prior information on the solutions. For instance, physiological insight on the human anatomy leading to specific joint limits, combined with a constraint to maintain balance along additional information on the environment were used to sketch sitting poses [START_REF] Lin | A sketching interface for sitting-pose design[END_REF]. The poses can also be constrained to lie in the space of "natural poses", extracted from a database of human motion [START_REF] Wei | Intuitive interactive human character posing with millions of example poses[END_REF][START_REF] Choi | Retrieval and visualization of human motion data via stick figures[END_REF].

These assumptions could be used to resolve ambiguities within the line of action framework we provide. However, they would also make the tool less flexible. In this work, we propose the alternative approach of constraining skeletal transformations to lie in the viewing plane. This way, users are free to exaggerate bending angles and making unbalanced poses, adding drama to their figures-as often done by cartoonists. On the other hand, the user has to turn the camera in order to "edit depth"-while previous edits orthogonal to the viewing plane are kept unmodified.

The line of action is covered by many textbooks and tutorials on drawing such as [START_REF] Lee | How To Draw Comics The Marvel Way[END_REF][START_REF] Blair | Cartoon Animation[END_REF][START_REF] Hart | How to draw animation[END_REF][START_REF] Brooks | The Complete Cartooning Course[END_REF][START_REF] Doucet | Line of action-mickey's christmas carol[END_REF]. Although they explain the concept quite well, they never give a formal definition that can be translated into mathematical terms. The books recommend drawing a single line, and then populating the body around it. They strongly advise using only C and S-shaped curves-which tend to produce more aesthetically pleasing, and also more physiologically plausible poses. The line has multiple purposes, best summarized as: "You can think of it as the back bone of a character or just as the imaginary line that dictates how the body will move" [START_REF] Doucet | Line of action-mickey's christmas carol[END_REF]]. In other words, it can be used to render explicit the state of mind of a static character, but also to convey upcoming motion as clearly as possible. In this work, we give a formal definition of the LOA and formalize how it relates to the character's body. Following our observations of these textbooks, we use the line to drive both the positions and tangents of a sub-chain of bones within the character body.

To our knowledge, the first use of lines of action in computer animation was for post-processing animations; i.e. to automatically exaggerate motions by stretching the character's limbs along the line [START_REF] Noble | Automatic expressive deformations for stylizing motion[END_REF]. Similarly, [ Öztireli et al. 2013] consider lines of action for creating new poses and animations, but their focus is on deforming objects based on curving the bones. In contrast, our approach uses LOA to pose standard articulated models, focusing on the issues in making the LOA an effective interface. A popular tool to pose characters with lines is the IK spline [Autodesk 2009]. However, in contrast to our goals, IK splines are 3D curves, are controlled by position targets only, and their mapping to the bones has to be set manually. Lastly, a very intuitive interface based on gestural 2D strokes was used to drive character animations [START_REF] Thorne | Motion doodles: an interface for sketching character motion[END_REF]]. The method focuses on trajectory and timing specifications to select and combine pre-computed motion clips. In contrast, the lines we draw specify poses, and our tool can be used to create new animation clips.

Overview

This work introduces a direct and natural method for character posing inspired by the way humans tend to sketch shapes in early stages of design. The line of action (LOA), an aesthetic curved stroke, allows the user to quickly specify the character's global shape in a single hand gesture. The pose can then be refined from other viewpoints and the remaining body parts can be posed using so-called secondary lines, as suggested in [START_REF] Abling | Fashion Sketchbook[END_REF]]. The method is fully automatic; the user does not need to specify any correspondence between the LOA and the character's body, nor to manually pre-select bones.

Although sketching curves in screen space shares similarities with the 2D stick figures approach introduced by [START_REF] Davis | A sketching interface for articulated figure animation[END_REF]], using a single line of action as an abstraction of the skeletal parametrization introduces several problems of its own. First, an LOA is a smooth curve, restricted to a specific set of shapes, namely S and C shapes. Second, by looking at many examples such as those in Fig. 2 or Fig. 4, we notice that the LOA dictates both position and tangent targets, but for only a specific set of bones in the character. Based on these insights, we derive an optimization problem in Section 2 to pose 3D characters from 2D sketches. This model is sufficient to handle perfectly drawn lines of action whose lengths' are the same as the characters' in screen space. In many cases however, the correspondence is ambiguous; lines are drawn too small or too large locally or globally. A good example is the S-shaped curve often used for feminine stances shown in Fig. 4. We can see an area near the pelvis which seems to correspond to none of the bones in the body. We propose an automatic method to solve this problem. At each step we solve for an optimal reparametrization of the correspondence between the line of action and the bones-which we describe in Section 3. Also, the "3D pose from 2D sketch" problem is well known for being under constrained; not enough information is provided by the sketch to fully specify the solution. Since one our goals is to allow expressive, sometimes exaggerated and non-realistic poses, we make only minimal assumptions on the character's physiology. To constrain the problem properly while allowing flexibility in the pose, we propose to constrain the character transformations to lie in the viewing plane (Section 5). Consequently, the user has to rotate in order to "edit depth" without ambiguity. This allows the user to draw poses not restricted by physiological assumptions such as balance, joint constraints, or pose priors-although such constraints are also compatible with our framework. Finally, we show results of our tool in section 6 and accompanying video.

The Line of Action

The line of action's description often differs from one artist to the next. To formalize the concept, we provide a definition which, in our view, spans a wide range of cases. A first aspect is the range of possible shapes an LOA can take, and the second is its relation to the character's pose. Both are discussed below.

LOA-shaped curves:

In most textbooks on drawing, the LOA is restricted to the family of Cand S-shaped 2D curves-which we refer to as LOA-shaped curves. This restriction can be a useful for two reasons: it constrains the poses to aesthetically pleasing shapes, and to physiologically plausible poses. We define C and S-shaped curves as having zero and one inflexion points respectively. To naturally reduce the variability of the line, we fit a cubic Hermite curve, defined by two end points and two tangent vectors. This is done as follows: the end-points of the Hermite curve are positioned at the extremities of the stroke. We then minimize the squared distance between sample points along the stroke and the curve positions, w.r.t. the angles and lengths of the tangent vectors, represented in polar coordinates. This approximation is sufficient for representing the usual C and S-shaped curves. More freedom to the user is given by fitting more curves to the input stroke; in our case two curves are used, each capturing half of the user's stroke.

Using the LOA to drive the pose: We use the LOA to dictate the positions and tangents in screen space of a subset of the character's bones-which we call the body line. For instance, we often see the line setting the bones going from the head to a foot-the yellow bones in Fig. 3. Our definition of the body line is as follows:

The body line is defined as a maximal, connected linear chain in a character's skeletal kinematic tree.

Maximal: Most body lines start and end at extremes of the full kinematic tree. Restricting our definition to maximal lines leads to a set of 10 possible kinematic chains. Different ways of choosing between these body lines are discussed in Section 4. For the sake of clarity, we will continue our description supposing a classic body line, i.e. a chain between the head and one of the feet.

Connected: The bones are children or parents of one another in the full kinematic tree. For instance, if the LOA goes from one foot to the other, it also includes the two lower and upper legs.

Linear: Each node between two bones, that are not the extremities, are of degree 2. The extremity nodes are degree 1 nodes; they have only one parent or one child.

There are two important aspects to note about the body line in regard to the line of action. One is the importance of shape: it implies that tangents should prevail over positions. It seems more important for the character to match the shape of the LOA than having the character's body line be close to the line of action. This will be reflected by a tangent matching constraint in our model.

The second is that some bones seem to be more important, and to align more with the curve. For instance, in Fig. 4, the lower leg, the spine and the head seem to be matching the shape of the curve, while the upper leg, the pelvis and the neck account for less. Although we describe a general model, we find important to emphasize that these principal bones will weigh more in practice. Both On the right is the result after solving the LOA problem (1). In the top image, the body line goes from the head to a foot. And in the second, it goes from a hand to a foot.

the sketched LOA (2D) and the character's body line (3D), can be seen as one dimensional parametric curves. For the sake of clarity, we use the same parameter s for both curves. In Section 3, we will make our solution more general by finding an optimal warping function between bone coordinates and curve coordinates.

Let the position in screen space of the line of action be x loa (s), the bones' positions x b (s) and their position in screen space Pvpx b (s), where Pvp is a view and perspective projection transformation in homogenous coordinates.

Following our definition, the best pose x b (s) of the character minimizes the shape difference between the two curves. Importantly, shape is defined here both in terms of positions and tangents. We thus solve for:

min x b (s) s Ex(s) + E T (s) ds, (1)
Ex(s) = λx(s) Pvpx b (s) -x loa (s) 2 , E T (s) = λ T loa (s) T b (s) -T loa (s) 2 ,
where T b (s) = ∂Pvpx b ∂s (s), T loa (s) = ∂x loa ∂s (s), and T denotes unit tangents. The terms λx(s) and λ T (s) are used to emphasize tangents over positions, and to emphasize the importance of specific bones over others, e.g. the principal bones. For instance, for the classic head-to-foot body line, the lower leg, the spine and the head usually weigh more than the upper leg, pelvis and neck. We minimize w.r.t. to the full body line x b (s), using a discrete version of the problem given in the appendix.

Automatic Correspondence

Driving the character's body line from a single stroke, requires identifying which part of the body corresponds to which part of the stroke. In Section (2), we supposed a known mapping, expressed through a common coordinate s, but this assumption does not hold in practice: the user-sketched stroke is often ambiguous and gives no information on which part of the LOA the bones should correspond to. Moreover, two simple ideas for setting this mapping just do not work:

• A first idea would be to match each projected bone to the nearest point on the LOA. Although this works well for a character that is already "near" the curve or already well aligned with the curve, it falls short for ambiguous Additionally, the closest point can lead to unexpected character configurations as shown in Fig. 4; we see the upper body sliding to the middle of the line of action when we would expect it to remain in the upper part of the LOA.

• A second idea is to consider the length of the bones. However, the sketched line is often drawn imperfect with local segments having a different length than the projected body line's length.

We even found cases where some parts of the LOA correspond to nothing in the body, as with the S-shaped curve shown in Fig. 4-with a jump visible near the pelvis.

We address these problems by finding an optimal spatial warping between the body line coordinate s and the line of action coordinate w(s). We formulate our approach on the closest point Ex(s, w(s)) objective, but add regularization on the warping function w(s) based on the assumptions of rigidity and connectivity of the skeleton's tree structure. Note that rigidity here is used to find the mapping, but the bones can be made to bend in a final step simply by breaking them into smaller ones.

In short, assuming the bones are rigid implies they should avoid covering areas on the LOA with high curvature which translates into constraint Eκ(w(s)). To allow jumps in the solution (see Fig. 4), the warping function is defined piecewise, each corresponding to a rigid bone interval on the body line denoted Ai. Assuming a skeletal tree structure, the associated warping functions should preserve the parent-child relations established by the skeleton; bones should not overlap or go too far from one anothermaking a connectivity constraint EC (w(s)) relevant. The final optimization problem is as follows:

min w(s) s
Ex(s, w(s)) + Eκ(w(s)) + EC (w(s)) ds , (2)

Eκ(w(s)) = λκ ∂ 2 x loa (w(s)) ∂ 2 s 2 , EC (w(s)) = λC χi+1(s)w(s) -χi(s)w(s) 2 ,
where the term Ex(s, w(s)) is the same as in problem (1), but here the line of action coordinate is warped w(s). The term χi(s) is a step function worth 1 at end points of the rigid bone interval Ai and 0 elsewhere. Note that we optimize for the warping function on the fly; i.e. at each step when solving problem (1), we solve problem (2) for an optimal warping function. Between successive steps, convergence for problem (2) is quite fast when initialized with the previous step solution. More details on our discretization are given in the appendix.

Selecting the Body Line

The line of action can be used to modify different sets of bones, or body lines. Sometimes the body line includes an arm, sometimes both legs, while in other cases it includes only the arms-as with secondary lines [START_REF] Abling | Fashion Sketchbook[END_REF]].

Our definition of the body line was designed to include all of these cases. Due to the maximality criterion, we only have 10 body lines for a humanoid: head to left hand, head to left foot, head to right foot, head to right hand, left hand to left foot, left hand to right foot, left hand to right hand, left foot to right foot, right hand to left foot and right hand to right foot. Selecting the right body line in the right viewpoint can be viewed as an artistic choice. Because we wanted to offer a simple interface, we propose an automatic way of selecting the "most appropriate" body line for a given line of action and viewpoint.

We propose two interfaces from which the user can choose in order to select a body line. The first is manual. The user either draws a first line close to the bones he wishes to include, or he selects the body line from a tab. Finally he draws the line of action. With the second interface, the user draws the target line of action directly and the system automatically chooses the body line with the smallest energy (1) from the 10 possible candidates. When two body lines have similar energies, we take the closest one in view space.

Resolving Depth Ambiguities

The line of action problem (1), which seeks a 3D pose from a 2D sketch given a viewpoint, is under constrained; many solutions (3D poses) exist for the same 2D line of action. For instance, we can simply imagine a line drawn shorter than the character in screen space. It could mean the character has translated away from the camera, but it could also mean that he has crouched. Since we don't make any assumptions other than the piecewise-rigidity of the skeletal structure-in order to allow non-realistic poses-we need additional constraints to resolve depth ambiguities.

Using the 2D LOA to drive the character in all three dimensions can be confusing for the user due to the unexpected character configurations; several 3D poses can be a minimum to problem (1). To "remove" depth, we constrain the character's transformations to lie in the viewing plane-which also reduces ambiguities when solving problem (1). We parametrize each bone rotation to a single axis-angle component θi as illustrated in Fig. 5. The corresponding axis is the camera's viewing direction projected onto the floor plane.

The translations are also parametrized to lie in the viewing plane, i.e. along the 2-components u, v of the viewing plane. Finally, we solve the LOA problem (1) with respect to θi and u, v.

With our approach, the user can separately sketch and rotate the camera to "edit in depth". As an example, if he wants to bend the character's knees, he has to look at the character sideways and draw a curve indicating a bent leg as shown in Fig. 5. Working only in the viewing plane removes a lot of the ambiguities-both for solving the problem and for the user as well. Without this constraint, we would necessarily have to supplement the problem with other constraints such as ensuring poses lie in the space of prior "natural poses" [START_REF] Wei | Intuitive interactive human character posing with millions of example poses[END_REF]. Although it would be very useful to do so, especially for beginners, it would also restrict the space of poses that can be created. The system used in this paper gives entire freedom to create non-realistic and exaggerated poses with arbitrary sequences of sketching and rotation gestures.

Results & Discussion

As a proof of concept, we evaluated our method with the task of reproducing static poses from photographs, cartoon images, and finally, example motions from video frames. The poses were then exported to Maya for keyframe interpolation. The implementation is a stand-alone software that offers all the necessary tools for drawing multiple lines of action in multiple viewpoints, and computing the corresponding poses at interactive rates. In our implementation, we solve the discrete version of problem (1)-discussed in the appendix-using gradient-based local optimization. At the start of each step, we solve the discrete version of problem (2) for an optimal reparametrization. In both cases we optimize until the gradient's length is smaller than a small criterion. The optimization is interactive and can be visualized on-line.

All the poses and animations were made by the authors of this paper-who are not animation experts-using lines of action only.

For the evaluation, we measure the total number of user gestures (strokes and rotations), as well as the actual time used to generate a satisfactory result. We start from a neutral pose, then draw all the strokes necessary to complete the pose or sequence of poses.

Then we asked a professional animator (with five years experience) to reproduce the same poses using IK widgets in Maya. We measured the time required to finish the same task and summarized our evaluations in a table. Our preliminary evaluation demonstrates that all the example poses could be generated within reasonable times. Figure 6: Number of strokes, rotations and time taken using lines of action v.s. using 3D IK widgets in Maya.

While we have not performed a complete comparison, we extrapolate that our tool can be used to produce rough poses and keyframes in a fraction of the time it would take with existing tools. The benefits seem especially spectacular when considering longer animation sequences like the Muybridge sequence-shown in Fig. 9 and the accompanying video. It took a total of 6 min. to make using lines of action while it took 22 minutes to a professional artist using traditional 3D IK widgets in Maya.

One of the reasons for this speedup is the reparametrization described in Section 3. Thanks to this adjustment, we can now sketch bent knees and elbows more freely without having to produce an accurate stroke each time. This is beneficial in reducing the number of strokes required to produce the poses-not to mention poses like the feminine stance (Fig. 1) that are nearly impossible to reproduce without a reparametrization. But the foremost benefit of our method is the fact that beginners could actually create expressive poses using a line of action, while there is no way in Maya for beginners to easily create poses that exhibit an aesthetic curved shape.

The line of action is particularly useful to create dramatic poses that convey full body expressions more clearly.

In practice we found that most keyframes were made by sketching 3 strokes: one to specify each half of the full body, then a third to adjust the arms using a secondary line. Moreover, it is worth noting that once a pose is created, variations of the same pose can be easily tried-out by drawing slight variations of the initial line. This explains why making a basic walking pose (Fig. 10) takes only 20 seconds while making a dramatic punch destined for comics like the one in Fig. 11 takes minutes; we tend to redraw the same stroke over previous ones to perfect the pose and make it more dramatic.

The poses and animations shown in this paper and video were made using lines of action only. However, there are many poses that just cannot be fully made using lines of action only-insofar the line is restricted to S-and C-shaped curves (see Section 2). Sharp corners such as completely bent knees cannot be reproduced using a limited number of smooth parametric curves. On the other hand, this option to constraint the user to make more consistent poses at the coarse level-a feature useful for beginners. Of course, LOA constraints can be combined with other techniques better suited for small scale detail editing.

A limitation of our viewing plane constraint is that bones can no longer twist when orthogonal to the viewing direction. We explored the possibility of letting the bones rotate around their own direction axis, but have not yet established an intuitive way of controlling this feature. Up to now we used the secondary lines for the arms or legs in order to twist the upper and lower bodies. For instance, this was used to produce Fig. 11. Future work could investigate considering the length of each arm or leg in order to estimate the amount of twist required.

Conclusion & Future Work

We introduced an intuitive interface for expressive character posing based on the line of action metaphor. We gave a formal definition of the line and provided solutions to several ambiguities related to this abstraction of the body; namely correspondence and depth ambiguities. Our method enable the user to rapidly create expressive poses ranging from elegant fashion stances to exaggerated cartoon actions. We briefly compared our method with traditional authoring software and found it was significantly faster for coarse pose design-especially for expressive poses that are challenging to create using traditional animation techniques. This new tool for character posing seems promising and opens several questions for future investigation.

One of the advantages of the optimization problem (1) is that we can easilly add additional constraints. For instance, when the line intersects the floor, we can handle a foot or hand placement IK target. This simply corresponds to an additional constraint in problem (1). We extended the line of action concept to several kinds of lines, but we restricted our definition to maximal chains (Section 2). It would be convenient to allow individual body parts like a leg or an arm to be sketched individually to refine poses without having to sketch a maximal line each time. On the other hand, our definition of the line is quite general and we could apply it directly to other morphologies like quadrupeds. An LOA could be formed from the head going through the spine and finishing at a foot. Secondary lines would go from one leg to the other. Also, we could think of refining small scale details like hands following the same logic. A line could be used to set one of the fingers to the thumb, or one finger to the elbow.

Our goal was to make the line of action as simple as possible for the user. For this reason we proposed fully automatic solutions to the problems we faced; namely the correspondence issue (Section 3) and the body line choice (Section 4). However there are cases where the user might want a different correspondence, or a different body line to modify. A semi-automatic approach could be envisioned here were the user would scribble ( [START_REF] Noris | Smart scribbles for sketch segmentation[END_REF]) or sketch small primitives to help both resolve the correspondence and specify the choice of internal body line. These insights are left for future work, but seem compatible with our framework; they could be inserted as additional constraints in problems ( 2) and (1) respectively. The LOA curve position is parametrized as a cubic Hermite curve (see Section 2) and can be sampled between nodes as: x loa (wi,j), where each node has a warped coordinate si → wi, and wi,j = wi + j n are the coordinates between the nodes. By discretizing the partial derivatives in problem (1) in finite difference, the discrete problem can be written as:

min Q,xr m+1 i=1 E x i + E T i E x i = λ x i n + 1 n j=0
Pvpx b (si,j)x loa (wi,j)

2 E T i = λ T i n + 1 n j=0
T b (si,j) -T loa (wi,j) 2 T b (si,j) = Pvpx b (si,j+1) -Pvpx b (si,j) ∆sj

T loa (wi,j) = x loa (wi,j+1)x loa (wi,j) ∆sj

where ∆sj = 1 n ∆si are subintervals' lengths. The matrix Pvp is the view and perspective projection transformation in homogenous coordinates and T = T T denotes unit tangents. The principal bones can be weighted more through λ x i and λ T i with the later being larger than the former to emphasize tangent matching over position matching. Note that in our implementation, the rotations are parametrized to a single axis-angle component θi constraining the transformations to the viewing plane and reducing the number of degrees of freedom to solve for (see Section 5). Finally, we minimize w.r.t. Θ = {θ1, ..., θm} and xr.

The warping function is piecewise-defined-one segment for each bone in the body line-accounting for potential jumps in the solution (see Section 3). We discretize these segments into two components W = [w 0 1 , w 1 1 ], ..., [w 0 m , w 1 m ] , one for each extremity of the bone. This means nodes with two edges also have two coordinates; we average both to sample the curve at a node i as wi = 0.5(w 0 i + w 1 i-1 ). By discretizing the partial derivatives in problem 2, a discrete version of the problem is:

min W m+1 i=1 E x i + E κ i + E C i E x i = λ x n + 1 n j=0
Pvpx b (si,j)x loa (wi,j)

2 E κ i = λ κ n + 1 n j=0
x loa (wi,j+1) -2x loa (wi,j) + x loa (wi,j-1) ∆sj 2 2

E C i = λ C w 0 i+1 -w 1 i 2
where λ x is to chose the closest point on the line while the terms λ κ and λ C are used to emphasize rigidity of the bones and their connectivity. Finally we minimize w.r.t. the warping function W on the fly, between steps for problem (1).

Figure 2 :

 2 Figure 2: Our definition of the line of action:The line is a target for a subset of the kinematic chain from a particular viewpoint. The bones are generally a lower leg, the spine, the head and sometimes the arms. ©The Estate of Preston Blair,[START_REF] Blair | Cartoon Animation[END_REF]].

Figure 3 :

 3 Figure3: In yellow, on the left, are the bones forming the body line that will be deformed according to the sketched line of action. On the right is the result after solving the LOA problem (1). In the top image, the body line goes from the head to a foot. And in the second, it goes from a hand to a foot.

Figure 4 :

 4 Figure 4: Left: A classic S-shaped curve used for feminine stances. It is difficult to know which part of the line, the bones should map to-not to mention parts near the pelvis that correspond to none of the bones. The three last images are: the initial pose, the result parameterization using the closest point only and finally, our result with curvature-based regularization. Left figure: ©Ben Jelter.

Figure 5 :

 5 Figure 5: Bone rotations are constrained to the viewing plane. They are parametrized as a single axis-angle component. The axis is the camera's viewing direction projected onto the floor plane. From left to right is the initial pose, a side view with the line of action stroke, the result after optimization, and a frontal view of the final pose.

Figure

  Figure Num. Strokes Time Time & Rotations LOA Maya Walk (Fig.10) 4 st. 1 r. 20 s. 90 s. S shape (Fig.4) 3 st. 0 r. 20 s. 120 s. Hero Punched (Fig.11) 6 st. 2 r. 90 s. 150 s. Hero Punch (Fig.11) 7 st. 4 r. 2 m. 3 m. 30 s. Animation Dancers (Fig.8) 6 st. 1 r. 30 s. 2 m. 45 s. Cartoon swing (Fig.7) 7 st. 1 r. 45 s. 2 m. 30 s. Muybridge (Fig.9) 34 st. 1 r. 6 m. 22 m.

Figure 11 :

 11 Figure 11: Inspired by comics, these characters are posed fighting. In each of these poses, we can see the overall body expression controlled by a line of action.

Figure 7 :

 7 Figure7: The main poses of a cartoon character swinging a bat take only 8 strokes and 1 rotation. The overall process took less than a minute to finish compared to several minutes with traditional 3D IK widgets.

Figure 8 :

 8 Figure 8: Dance poses are often LOA-friendly with the body exhibiting an aesthetically curved shape.

Figure 9 :

 9 Figure9: A complete Muybridge sequence can be reproduced in 6 minutes by drawing lines of action in a single view. We show, for each frame, the last line of action that was drawn.

Figure 10 :

 10 Figure 10: Starting from a neutral pose, only 4 strokes are required to get started with a walking animation. Note the so-called secondary lines are used here to place the arms.

Acknowledgements

We thank Laura Paiardini for help with Maya scripts, rendering and animation comparisons. We also thank the anonymous reviewers and Mike Gleicher for their useful comments and suggestions. This work was partially funded by the ERC advanced grant EXPRES-SIVE.

Appendix

We parametrize the bone position function x b (s) as a chain of piecewise-rigid bones: a single root position at the pelvis xr with a quaternion rotation Q = {q1, ..., qm} for each bone. The chain holds m + 1 nodes and we denote x b i the bone position at node i on the body line. The body line coordinate s is discretized into intervals-one for each bone-starting at si and covering ∆si = si+1si space. We denote si,j = si + j n a coordinate between two nodes and build a piecewise continuous bone position function by linearly interpolating between the nodes:

, for a subinterval of n + 1 samples (n = 4).


Mickaël Duruisseau

Romain Rouvoy

Automatic Inference of Roadmaps from Raw Mobility Traces

Keywords: clustering, roadmap inference, smartphone, urban sensing clustering, inférence, carte routière, smartphones

The popularization of smartphones in daily life offers numerous opportunities in terms of urban sensing. More and more users are ready to share certain information as part of scientific research, including their GPS location. From these mobility traces, we developed a roadmap inference algorithm using raw mobile data supplied by users of smartphones. This algorithm can generate a map composed of oriented routes, which are annotated by a certain amount of metadata.

Introduction

Nowadays, there are many maps produced manually or via the use of special cars, such as Google Maps1 or OpenStreetMap 2 . Although these maps are pretty accurate, they are not always upto-date, and requires continuous interventions to update them.

Inferred maps, however, may be used to detect changes to the road network or errors in existing maps, as well as being updated in real-time. Map inference may also be used to produce custom maps for certain classes of travelers, such as pedestrians or cyclists, by feeding it with data from different sources.

However, inference requires a lot of data to be effective, which was a limitation until recently. The GPS3 traces needed for producing such maps are becoming increasingly available due to the popularity of smartphones and GPS navigators, which capture vast volumes of user location traces on a daily basis.

The main goal and motivation of this work is to be able to infer roadmaps on a large scale, which is not possible when using pictures [START_REF] Biagioni | Map inference in the face of noise and disparity[END_REF][START_REF] Biagioni | Easytracker: automatic transit tracking, mapping, and arrival time prediction using smartphones[END_REF] Also, keeping the road bearing or accuracy is not easy with images, nor is the ability to extract several lanes on the same roads.

Additionnally, our goal is to update this map in real-time, just by adding data, collected via crowd-sensing systems, for example. This can be particularly useful when a road is closed due to roadworks, hazard or anything that could block it. Moreover, depending of the data, we can infer pedestrian ways, bike trails or anything like that.

The contributions described in this paper are:

1. A three dimensional clustering algorithm for GIS 4 datasets, based on coordinates and orientation, 2. A weighted distance calculation using distance, bearing and accuracy, 3. A graph database storage to represent roads by relationships between nodes, 4. A road-maps inference algorithm and several filtering algorithms based on this database.

State of the Art

In this section, we introduce various roadmap inference approaches, which have been already reported in the literature.

K-Means Clustering

This family of algorithms, originally described in Edelkamp and Schrödl [START_REF] Edelkamp | Computer Science in Perspective[END_REF], consists in distributing a serie of cluster seeds at different locations based upon input trace data. These seeds are placed following constraints on bearing and distance and are used as initial guesses for a standard K-Means clustering algorithm [START_REF] Edelkamp | Computer Science in Perspective[END_REF]. Once the clusters have been built, they are linked to form road segments based on the patterns that the raw traces used to follow initially.

Edelkamp and Schrödl [START_REF] Edelkamp | Computer Science in Perspective[END_REF] The paper by Edelkamp and Schrödl, is the original k-means based method, but it also refines the road network model by fine-tuning the location of intersections, representing the road center-line by a fitted spline rather than a series of straight lines.

Schroedl et al. [START_REF] Schroedl | Mining GPS Traces for Map Refinement[END_REF] This paper describes a process to refine the intersection geometry and model individual lanes. This process is accomplished by identifying intersections and their bounding boxes, and group traces by entry and exit points. Then, a spline-fitting technique produce the final turn-lane intersection.

Agamennoni et al. [START_REF] Agamennoni | Robust inference of principal road paths for intelligent transportation systems[END_REF] This paper uses an approach similar to Schroedl et al. [START_REF] Schroedl | Mining GPS Traces for Map Refinement[END_REF], but focuses to extract principal road paths.

Trace Merging

Trace merging algorithms iterate over each GPS trace recorded by adding them to a map. For each segment from the trace, a test is performed to determine if it already exists on the map or if it is "similar enough". In this case, the segment is merged into the existing one, otherwise it is added as a new segment. When a segment is added to an existing one its weight is increased.

Segments with a weight below a certain threshold are removed.

Niehoefer et al. [START_REF] Niehöfer | GPS Community Map Generation for Enhanced Routing Methods Based on Trace-Collection by Mobile Phones[END_REF] The method described by Niehoefer et al. [START_REF] Niehöfer | GPS Community Map Generation for Enhanced Routing Methods Based on Trace-Collection by Mobile Phones[END_REF] evolves the merging approach by adjusting the position of road segment when merging a new trace. This technique allows the location of road segment to be refined as more traces are added. This paper also describes a way to automatically classify road types, such as streets or highways.

Kernel Density Estimation

Kernel Density Estimation (KDE) algorithms begin by discretizing the geometric space covered by the set of GPS points into a grid of pixels, recording the approximate density of samples for each cell. A threshold is applied to the density map in order to produce a binary representation of the road network, and the center-line is extracted using a variety of image-processing methods.

Biagioni et al. [START_REF] Biagioni | Map inference in the face of noise and disparity[END_REF] Kernel Density Estimation is heavily dependent on the threshold parameter used to extract the boolean map. This introduces the tradeoff between map accuracy and road coverage. Since it is such a sensible value, it is proposed by Biagioni and Eriksson [START_REF] Biagioni | Map inference in the face of noise and disparity[END_REF] that no single threshold can produced the desired map thereby a methodology for constructing a map with various thresholds is proposed.

Contributions

Input Datasets

In this work, we used two different datasets. The first one is provided by the state of the art, and we use it to compare the accuracy of our algorithms. The second is a dataset we created with the help of a colleague who works on smartphone energy consumption.

Chicago Bus Shuttle This dataset is used by Biagioni et al. [START_REF] Biagioni | Map inference in the face of noise and disparity[END_REF][START_REF] Biagioni | Easytracker: automatic transit tracking, mapping, and arrival time prediction using smartphones[END_REF] in their comparison and evaluation of the state of the art. Almost one 120, 000 locations are reported, but with only their locations and timestamps. Furthermore, each location (except the first and the last ones) is linked to the next and previous locations in the mobility trace.

Villeneuve d'Ascq Daily Trips As described earlier, this dataset was created by Inria and contains more than 22, 000 locations. Compared to the Chicago one, this dataset provides more metadata for every location, especially the accuracy and the bearing. This information is very useful to filter spurious locations before applying our algorithm.

Clustering Algorithm

The first contribution reported in this paper is a clustering algorithm used to group raw locations into groups of correlated locations. We first describe the data structures used to store and manipulate the location traces we extracted from the datasets and then report on the algorithms defined to process these structures.

Data Structures

We use several structures to represent data. In our program, we employ Scala5 as main language, and choose to implement those structures as classes:

GPSNode -A raw location, it contains all the information retrieved by a GPS fix, plus the trip and cluster to which the node belongs, and its neighbours nodes.

GPSCluster -A group of correlated nodes, its own location (latitude, longitude and bearing) is updated whenever a GPSNode is added and computed as the geometric center of the nodes it encloses.

Trip -A list of consecutive nodes, linked to each other.

Figure 1 is the UML6 class diagram of our project. The first step of our procedure is to create trips (Trip) filled with nodes (GPSNode). We parse the dataset files, create a GPSNode per line with all the data provided, and make links between them, in order that nodes are linked by their properties next and prev. For example, N 1 .next = N 2 , N 2 .prev = N 1 , as depicted in Figure 2.

Those links are used to loop over a trip or a road, and to calculate the velocity of a node, together with its bearing, if it was not already set (depending of the dataset). Then, we use Algorithm 1 to create the location clusters from the raw locations, this algorithm iterates over the list of trips and the locations each trip contains.

Additionally, we feed a database with all the locations, clusters and roads created by this algorithm.

Then, we apply several filters on the model to clean the clusters and inferred roads, as lated explained in Section 3.5.

Distance Computation

Our clustering algorithm uses a specific distance computation (function distance in Aglorithm 1), to create clusters depending on three parameters:

• The spatial distance, • The bearing difference,

• The location accuracy.

The distance formula we defined is customizable with a weight p. This weight allows us to change easily the importance of the distance or the bearing, depending of what we think will make better clusters.

distance(N 1 , N 2 ) = p × f d (d N1,N2 ) + (1-p) × f b (bd N1,N2 ) (1)
3D distance formula

In this formula, f d and f b are Gaussian functions parameterized by a threshold, they return 0.5 for a given value.

The result of that computation is a value normalized between 0 and 1. 1 meaning that the two input locations are exactly similar, 0 meaning that they are too far or with a high bearing difference.

In this paper, we consider that two nodes are close only if the distance is more than 0.5.

Spatial Distance

We use the Haversine Formula [START_REF] Sinnott | Virtues of the Haversine[END_REF] to compute the distance sd between two points. This formula gives great-circle distances between two points on a sphere from their longitudes and latitudes.

R = 6371.0//km dLat = rad(lat x2 -lat x1 ) dLon = rad(lon x2 -lon x1 ) lat1 = rad(lat x1 ) lat2 = rad(lat x2 ) a = (sin( dLat 2 )) 2 + (sin( dLon 2 )) 2 + cos(lat1) × cos(lat2) c = 2 * atan2( √ a, √ 1 -a) d = R × c sd = d * 1000//m (2)
The Harvesine formula Even if this formula uses an approximate value for the Earth radius, the results given are accurate enough for short distances. Once we have computed that distance, we apply an accuracy threshold, explained in Section 3.3.3. Then, we use a Gaussian function on the result, to normalize the distance into a value between 0 and 1, that represents the similarity between two nodes. RR n°8585 s = 5.0

f d (d) = e -(λ d ×d) 2 , λ d = log 2 s 2 (3)

Gaussian function for distance normalization

This formula is set by the value s, meaning that it will return 0.5 if d == s, more for a lower s, and less for a higher s.

In our algorithm, we set s to 5.0 because we wanted to be able to separate multiple lanes of the same road, for example on a highway.

It also means that our clusters will be constituted of the closest nodes with a maximum distance of 5 meters.

Bearing Difference

In addition to the spatial distance, we also use the bearing in our distance computation. It allows us to create oriented clusters and thus oriented roads. We apply the same procedure as for the spatial distance-i.e., we calculate the difference between the two bearings-and then apply a Gaussian function.

First, we have to retrieve the bearing of a node. There are two ways: either the bearing is returned by the satellite when it fixes the location, or we calculate it with the following formula.

Bearing between two locations Then, we obtain an angle in degrees between the two bearings, we want to get the difference between two angles, modulo 360. To do that, we use the following formula, where a and b are two angles in degrees.

To finish, we apply a Gaussian function, formula 6, as we did for the spatial distance, that returns the similarity between two angles. This function will return 0.5 for a given value s, more for a lower s, and less for a higher s. In our algorithm, we set s to 15.0, it allows us to detect any kind of change in the road direction.

The result of this function is then added to the spatial distance with a weight (the complement of p, (1-p)) to get the final result of our distance computation.

bd x1,x2 = 360 + b x1 -b x2 (mod 360) if (bd x1,x2 ) > 360 2 ) bd x1,x2 = 360 -bd x1,x2 (5)
Difference between two angles, modulo 360

f b (b) = e -(λ b ×b) 2 , λ b = log 2 s 2 s = 15.0 (6)
Gaussian function for bearing

Accuracy Threshold

Depending on several physical conditions, a GPS fix can be more or less accurate. For example, in large cities, urban canyons [START_REF] Schuessler | Processing raw data from global positioning systems without additional information[END_REF] can reduce the vibility of available satellites, and so the accuracy. Every fix is therefore supposed to be annotated with an accuracy.

We chose to be optimistic in our approach: we deduce the average accuracy of the two nodes to their distance, as shown in the formula 7.

d x1,x2 = max(0, sd x1,x2 - acc x1 + acc x2 2 ) (7)

Distance between two nodes, with accuracy

This method allows to merge nodes that are far from each other but also inaccurate, as they could be close if the fix was more accurate. This also helps to create cleaner roads, since we merge nodes and move the clusters according to them.

Database Storage

The database storage was one of the most sensible point to take care of, because we wanted our solution to be scalable, fast, and efficient. Thanks to our program architecture, we can easily change the database type we want to use, this allowed us to try our algorithms on several technologies.

At first, we tried to use a Scala R-Tree implementation, named Archery [START_REF] Archery | [END_REF], and to store our dataset in memory. It worked fine, but there was no data persistence, and the associated query engine was not powerful enough.

Then, we used some relational database, like PostgreSQL [START_REF]PostgreSQL[END_REF], with a spatial plug-in to easily store and request on geo data. This kind of database works well, but relationships between nodes are not easy to store neither to query.

Finally, we switched to graph databases, that fits all our requirements, as it can store nodes and edges, and have a fast and powerful query engine. We will see in section 3.4.1 more about graph database.

Graph Database

The main advantages of a graph database-compared to a relational one-is the ability to link nodes to each others, and be able to execute powerful query using those edges. For our research, we study several graph databases, we wanted one that was fast, scalable and with spatial support.

The first one we tried was Neo4J [START_REF]Neo4J -A graph database[END_REF], which comes with a very powerful query engine, and a data visualization website. Our algorithm was pretty fast with it, but despite those advantages, Neo4J does not handle very well concurrency, and its REST API is not fast enough.

Then, we used Mahout and Hadoop, but even if, at first, our algorithm was using MapReduce [START_REF] Robson | Clustering very large multi-dimensional datasets with MapReduce[END_REF][START_REF] Kim | DBCURE-MR: An efficient density-based clustering algorithm for large data using MapReduce[END_REF][START_REF] Fries | Projected Clustering for Huge Data Sets in MapReduce[END_REF][START_REF]Applying Hadoop's MapReduce Framework on Clustering the GPS Signals through Cloud Computing 1[END_REF], they did not fit our requirements. GraphX and Spark, however, was a good alternative. Fast and with a great scalability and concurrency, despite some bugs, the main problem was the limited API, which did not allow us to store multiple type of data, neither make fast requests or clusters.

Finally, we chose to use Titan DB [START_REF]RR n°8585 RESEARCH CENTRE LILLE -NORD EUROPE Parc[END_REF], that we will present in section 3.4.3.

Graph Structure

Beside the database, our graph structure did not change. We store nodes and clusters with all their properties, and relations between them.

There are three types of relations, as shown in Figure 3.

road -Represents a road segment, between two clusters inside -Means that a node is inside that cluster next -Links the nodes to their neighbours

C 1 C 2 C 3 N 1 N 2 N 4

Titan Distributed Graph Database

Titan [START_REF]RR n°8585 RESEARCH CENTRE LILLE -NORD EUROPE Parc[END_REF] is a scalable graph database that supports transactions and so concurrency. It also supports several storage and index back-ends, as well as geographic data. So, everything our algorithm needs is provided by Titan.

RR n°8585

It comes with a query language, named Gremlin [START_REF]Gremlin -Graph Traversal Language[END_REF], a powerful Java API, and a graph server called Rexster [START_REF]Rexster -A graph server[END_REF], which is useful to visualize and query data. With Titan, we had to describe the graph first. As presented in Section 3.2.1 and in Figure 3, we store GPSNode and GPSCluster, with all their properties, plus a property type, which indicates whether the node is a node (GPSNode) or a cluster. With that, we defined the three possible relations: inside, next, and road.

There was no need of a prev relation, as in the structure, because Titan allows us to get edges that comes in and out of a node. So almost every nodes have two next relations, one that goes out, and one that goes in, which can be used as prev relation.

We also created several indexes on our data. The first one is on the type property, and allows us to quickly get every nodes of a given type. The second is for the geo property, to get nodes close to a given location. And the third is a combination of the two first, type and geo, to get nearby clusters as used in Algorithm 1.

Filtering Algorithms

Once our algorithm 1 is completed, we obtain a large number of clusters and roads. Some of these clusters are incorrect or should be merged with a neighbour, and this applies also for roads. In this section, we will present four ways to improve data results, without loosing information.

Inconsistent Nodes

Even with our accuracy threshold (cf. Section 3.3.3), sometimes a cluster is created with one or two nodes, especially if the dataset does not provide a decent accuracy-it is the case for the Chicago dataset, for example. In this case, we apply a threshold filter on clusters. For every clusters, we calculate the average number of nodes in nearby clusters, and we delete the cluster if it has less nodes than half of this average (we chose to delete those with less than avg 2 , but this is easily customizable). This method is useful to delete clusters created from inaccurate nodes, even if accuracy is not set, and so remove the noise from GPS, or potential errors due to urban canyon or something else that reduces the number of available satellites.

RR n°8585

However, with data whose accuracy is defined, it is less useful, because inconsistent clusters are less likely to exist.

Duplicated and Loop Roads

Our next clean-up filter focus on roads, more precisely, roads that already exist, or that come and go in the same cluster. Duplicate roads could be hard to find if we had to iterate on every roads twice, but instead of doing that, we chose to iterate over clusters. For each cluster, we check if roads are unique or not-i.e., if an outgoing cluster has already been found-and we delete the associated road.

One good way would have been to set a road as unique when we defined the graph, but it was not possible with a many to many relation. The second-road loops-is easy to find and only requires one loop. For every road, we just have to check if the start node is the same as the final node. In this case, we remove the road, since there cannot be roads which come and go in the same cluster. Figure 6 shows these cases. Here, road 1 and road 3 are automatically removed by our second filter.

Triangle & Trapeze Paths

As shown in Figure 7, triangle paths appear when two nodes have a next in different clusters. This can happen if the time between two nodes is longer than usual (e.g., because of traffic).

Trapeze paths appear in the same conditions, but with an intermediate in addition, C 3 in the graph.

Generally, roads created because of that are "shortcuts" between two nodes, and are not correct. It can happen with just one (triangle) or two (trapeze) intermediates, but also with more. We chose to detect and delete only triangle and trapeze, as this algorithm can be applied several times, the remaining shortcuts would be removed on the next execution. In order to do that, we had to iterate over clusters. For each cluster, we get every clusters at a distance of one "road" (i.e., clusters are direct neighbours), and we check if it is not also at a distance of two or three. In this case, we remove the edge with the distance of one which is the incorrect "shortcut".

In Figure 7, a triangle is visible with clusters C 1 , C 2 and C 4 , the shortcut is road 1,4 and will be removed by our algorithm. It is the same for trapeze, the shortcut is road 1,4 , as a path between C 1 and C 4 , already exists. In Figure 8, a triangle is visible on the top right of the image, and another at the bottom, but less visible. Our filter deletes the shortcut road, that goes across the trees.

Square Paths

Our last clean-up filter is about square paths. As shown in Figure 9, it is possible to find roads that form a "square". In this case, since we cannot know which clusters are incorrect and we chose to merge the two clusters. Because our fusion updates the cluster location, the new cluster is at the weighted average location between the two merged clusters, according to the enclosed nodes. This means that if one cluster has one hundred nodes and the other only two, its location will not change much. However, if the two clusters are about the same node number, the new location will be between them. For example, if C 2 and C 3 have the same nodes number, their fusion would create cluster C 2,3 .

Then, we recreate the roads. To do that, we iterate over every clusters, and for each neighbour at a distance of one, we check if they have a same neighbour in common. In this case, if the distance between those two is not too high, we merge them. In our datasets, this cleaning is useful when a road has a part with less available satellites, leading to less accurate nodes, and possibly two "roads" on the same path.

As we can see in Figure 10, due to a lack of precision, two lanes were created, and then, the fusion centers the road. The top cluster was slightly down when merged.

Discussion

As mentioned in Section 3.1, we used two different datasets to test our algorithms. Our dataset provides the nodes accuracy, which is very useful in our distance computation to infer clusters. An inaccurate node will be added to a cluster even if it is slightly too far. Since a cluster location changes when a node is added, this merging does not add inconsistency to the cluster. Figure 11 depicts the resulting roads. The Chicago dataset, however, was more difficult to infer, due to the lack of accuracy. We tried several methods, as described in sections 3.3.3 and 3.5. Obviously, results are less accurate but remain correct.

All our formulas and algorithms are customizable, and we tested multiple values for each of those parameters in our work.

The main parameters are:

• Gaussian distance threshold, fixed at 0.5 meters,

• Gaussian bearing threshold, fixed at 15 degrees,

• Spatial distance weight, fixed at 0.5.

These parameters greatly affect the final result, and may be different depending on the dataset.

Conclusion

In this paper, we presented an updated state of the art, as well as a new method for road-maps inference. This technique consists in creating oriented clusters, and linking them to form roads. Compared to state of the art methods, this one is scalable, and stores many information about the clusters and roads, which bring new opportunities.

Indeed, our algorithms are just a proof of concept of roadmaps inference on a large scale, without using image processing. In particular, we can still improve and optimize our models, especially the distance computation and the clustering itself.

Furthermore, some new features could be added, such as stop or red lights detection, but also roads cutting to create routes and navigation helps. Activity detection of the user could also help us to separate the different types of paths, such as bike trails, roads, or pedestrian paths.

Figure 1 :

 1 Figure 1: Class diagram

Figure 2 :

 2 Figure 2: Links between two nodes

  dLat = rad(lat x2 -lat x1 ) dLon = rad(lon x2 -lon x1 ) lat1 = rad(lat x1 ) lat2 = rad(lat x2 ) y = sin(dLon) × cos(lat2) x = cos(lat1) × sin(lat2) -sin(lat1) × cos(lat2) × cos(dLon) b = deg(atan2(y, x)) + 360 (mod 360)

Figure 3 :

 3 Figure 3: Graph DB structure

Figure 4 :

 4 Figure 4: Rexster visualization of a cluster and its nodes

Figure 5 :

 5 Figure 5: Clusters filter

Figure 6 :

 6 Figure 6: Loop and duplicated roads

Figure 7 :

 7 Figure 7: Triangle & Trapeze paths

Figure 8 :

 8 Figure 8: Triangle path -Before (left) & After (right)

Figure 9 :

 9 Figure 9: Square paths

Figure 10 :

 10 Figure 10: Square path -Before (left) & After (right)

Figure 11 :

 11 Figure 11: Inferred roads

RR n°8585

http://maps.google.com

http://openstreetmap.org

Global Positioning System: http://www.gps.gov

Geographic Information Systems: http://en.wikipedia.org/wiki/Geographic_information_system RR n°8585

http://www.scala-lang.org

Unified Modeling Language: http://www.uml.org/ RR n°8585