Yohannes A Alamerew
email: yohannes.alamerew@grenoble-inp.fr

Daniel Brissaud
email: daniel.brissaud@grenoble-inp.fr

Evaluation of Remanufacturing for Product Recovery: Multi-criteria Decision Tool for End-of-Life Selection Strategy

Keywords: Remanufacturing, End-of-life strategy, Product Recovery, Multi-criteria decision tool

Remanufacturing is an end-of-life product recovery strategy whereby used products are restored to the original equipment manufacturer (OEM) standard and receive a warranty at least equal to a newly manufactured product. Current end-of-life product recovery decisionmaking approaches are focused on technical and economic factors neglecting crucial areas such as customer demand and legislative pressure which are pertinent in the decision-making process, more importantly in the case of remanufacturing. Additionally, there is lack of a holistic approach that adopts an inclusive methodology for analysing multiple conflicting factors: ecological, technical, economic, business and social aspects simultaneously.

This paper presents our research works on evaluation of remanufacturing for product recovery at strategic level. The paper identifies key end-of-life decision-making factors that aid for assessing the feasibility of remanufacturing. A Product Recovery Multi-Criteria Decision Tool (PR-MCDT) is proposed to evaluate the viability of conducting remanufacturing form an integrated point of view i.e. by simultaneously taking into account technical, economic, environmental, business and societal aspects. The main benefits of this tool are, a) the decision maker has the opportunity to consider key factors such as customer demand, legislative pressure and new technologies in the decision-making process, b) the tool helps to evaluate different end-of-life options under several often-conflicting criteria which are not necessarily quantitatively defined, c) the decision-making tool takes into account the preferences of the decision-maker in evaluating the feasibility of end-of-life strategies. An example is presented and discussed to illustrate the applicability of the tool for selection of product retirement option. The PR-MCDT is used by senor/middle management level to provide an early stage feasibility analysis of adopting remanufacturing for a general product type.

Introduction

The global crisis in resource scarcity, population growth and climate change impacts are placing pressure to ditch the traditional "Make-Use-Dispose" economic model and adopt "make, use, return" as our collective mantra by joining the circular economy. The implementation of extended producer responsibility (EPR) in new governmental legislation, together with the growing environmental and economic concern, demands that original equipment manufacturers (OEMs) to take care of their products after they have been discarded by the consumer [START_REF] Hosseinzadeh | An Optimization Model for Reverse Logistics Network under Stochastic Environment Using Genetic Algorithm[END_REF]. Product recovery has become increasingly important towards transitioning to a circular economy. Product recovery management aims to close the loop throughout the product life cycle [START_REF] Krikke 2 | On a medium term product recovery and disposal strategy for durable assembly products[END_REF].

Product recovery management (PRM) is the management of all used and discarded products, components, and materials to recover as much of the economic and ecological value as possible thereby reducing the quantity of discarded waste [START_REF] Thierry | Strategic Issues in Product Recovery Management[END_REF]. End-of-life product recovery strategies include Remanufacture, Repair, Recondition, Cannibalization, Redesign, Refurbish and Recycle [START_REF] Thierry | Strategic Issues in Product Recovery Management[END_REF] [START_REF] Jawahir | Technological Elements of Circular Economy and the Principles of 6R-Based Closed-loop Material Flow in Sustainable Manufacturing[END_REF]. All these end-of-life options are distinct from one another and selecting the best suitable product recovery option should take several factors into consideration [START_REF] Kumar | Value flow characterization during product lifecycle to assist in recovery decisions[END_REF].

Remanufacturing is one of the most important and emergent product end-of-life strategy [START_REF] Stewart | Moving forward in Reverse: A review into strategic decision making in Reverse Logistics[END_REF] for increasing resource efficiency and realizing circular economy [START_REF] Umeda | Modeling competitive market of remanufactured products[END_REF]. Remanufacturing is defined as "a process of returning a used product to at least original equipment manufacturer (OEM) performance specification from the customers' perspective and giving the resultant product a warranty that is at least equal to that of a newly manufactured equivalent" [START_REF] Matsumoto | Remanufacturing[END_REF].

Several studies have focused on evaluating the feasibility of conducting remanufacturing at strategic level. Dumande [START_REF] Dunmade | PLETS model: A sustainability concept based approach to product endof-life management[END_REF] developed PLETS model, a product life cycle extension technique for product end of life management. [START_REF] Sharma | 56)Benchmarking of product recovery alternatives in reverse logistics[END_REF] and [START_REF] Ilgin | Use of MCDM techniques in environmentally conscious manufacturing and product recovery: State of the art[END_REF] uses multi-attribute decision making approach to select the best product recovery alternative in reverse logistics. Chen and Wu [START_REF] Chen | Neural Network Model for Product End-of-Life Strategies[END_REF] applied extension of the End of Life Design Advisor (ELDA) tool using a neutral network model, while Gehin [START_REF] Gehin | A tool to implement sustainable end-of-life strategies in the product development phase[END_REF] developed a tool called Repro 2 designed to evaluate product suitability for remanufacture based on product profiles. [START_REF] Meng | Multi-objective optimization decisionmaking of quality dependent product recovery for sustainability[END_REF] developed a multi-objective optimization model to identify optimal product recovery solution. Strategic evaluation of remanufacturing is critical to ensure strategic decisions, which then promote success of the company [START_REF] Diaz | International Journal of Production Research Evaluating strategic remanufacturing supply chain decisions Evaluating strategic remanufacturing supply chain decisions[END_REF].

Current end of life product recovery decision-making approaches are centred on economic and technical factors neglecting other equally influential aspects which are pertinent in the decision-making process such as market demand, social trends and legislative pressure. This is highly significant in the case of remanufacturing. Additionally, there is lack of a holistic approach that uses an inclusive methodology to assess and evaluate the viability of conducting remanufacturing form an integrated point of view i.e. by taking into account technical, economic, environmental, business and societal aspects simultaneously.

The aim of this paper is to identify EoL decision-making factors and incorporate them into a holistic methodology to evaluate the feasibility of conducting remanufacturing. The viability of remanufacturing is evaluated with the relevant technical, economic, environmental, business and social criteria. The research objective of this paper is therefore to answer the following questions: - Which key factors should be considered in the evaluation of remanufacturing with respect to the relevant technical, economic, environmental, business and social criteria?  How to evaluate the feasibility of conducting remanufacturing holistically by analysing the different types of factors: ecological, technical, economic, business and social factors simultaneously?

The rest of the paper is organized as follows. Section 2 presents the research methodology used to answer the research questions. In section 3 the multi-criteria decision making approach is discussed and key decision-making factors used to assess the feasibility of remanufacturing are presented. Subsequently in section 4 the application of the method to a case is discussed. Finally, conclusions are drawn by summarizing the main findings of the study.

Research Methodology

As introduced in section 1, the goal of this research is to identify key decision-making factors and integrate into the proposed holistic decision-making tool that evaluates the viability of remanufacturing for product recovery by taking in to account ecological, technical, business and societal aspects simultaneously.

A comprehensive literature review was undertaken to identify key decision-making factors that used to assess the viability of conducting remanufacturing. Firstly, an exhaustive list of factors was presented and then the decision-making factors were sorted into main categories by the authors. Afterwards, factors from each category were evaluated based on literature review and expertise from G-SCOP laboratory. Based on the analysis, key decisionmaking factors were identified in regard to technical, economic, business, environmental and societal aspect and the most important factors are incorporated into decision making criteria. Based on findings from literature and feedback from expertise, most important factors pertinent to consider in the decision-making process are accentuated.

The multi-criteria decision making approach has been chosen as methodology to evaluate the feasibility of conducting remanufacturing at strategic level. An iterative and multi-level procedure is used for selecting an appropriate multi-criteria decision making methodology. The decision-making approach advised economic, environmental and societal indicators which will be integrated into the evaluation process. To show the application of the decisionmaking approach an illustrative example is presented based on a case study on an automotive engine. To simplify the complexity of the problem only the main components of the engine are considered.

Result and Discussion

Multi-Criteria Decision Tool

A Product Recovery Multi-Criteria Decision Tool (PR-MCDT) is proposed to evaluate the viability of adopting remanufacturing for product recovery. The six basic steps that grid the approach are as follows: (1) selection of potential end-of-life strategies, (2) scoping of end-of-life strategies, (3) selection of relevant indicators, (4) assessment of end-of-life strategies, (5) analysis and evaluation of EoL strategies, (6) refinement of strategies and final evaluation. Figure 1 shows the main steps of multi-criteria decision tool. MCDT is capable to consider product EoL selection from an integrated point of view i.e. by simultaneously taking into account environmental, technical, economic, societal and business criteria [START_REF] Kiritsis | Multi-criteria decision aid for product end of life options selection[END_REF]. The main benefit of this methodology comprises, the decision maker has the opportunity to consider constraints such legislation, new technologies and market demand. The decision-making approach also takes into account the preferences of the user in the evaluation process of end-of-life strategies. A brief description of each step of the tool is presented below.

i. Selection of Potential End-of-life Strategies

The definition of product recovery EoL strategies constitute the description of the product and associated potential EoL options. In this first step of MCDT approach, the decisionmaker identifies potential EoL product recovery strategies and is unlimited by any constraints.

Product recovery EoL options includes Repair, Recondition, Remanufacture, Cannibalization, Refurbish and Recycle. Table 1, presents a summary of main product recovery strategies. An EoL option is considered as a product recovery strategy, if it fulfils three main criteria: collection of used products, reprocessing of a recovered product and redistribution of the processed product [START_REF] Thierry | Strategic Issues in Product Recovery Management[END_REF].

A potential product recovery EoL strategy is a possible candidate for evaluation and comparison during the decision-making process [START_REF] Bufardi | Multicriteria decision-aid approach for product end-of-life alternative selection[END_REF]. In multicriteria decision literatures, list of potential candidate strategies are generally called alternatives or actions [START_REF] Roy | Multicriteria Methodology for Decision Aiding[END_REF]. A functional description of the product is decisive for the remanufacturing company to be able to achieve high level EoL treatment. The description of the product provides relevant information regarding the characteristics of the product as well as its functional use for the consumer [START_REF] Lamvik | The AEOLOS methodology[END_REF].

Definition of Product Recovery Strategies

Source Remanufacture is an end-of-life product recovery strategy whereby used products are restored to the original equipment manufacturer (OEM) standard and receive a warranty at least equal to a newly manufactured product.

[20][21] [START_REF] Rose | Design for Environment : A method for formulating end-of-life strategies[END_REF] Recondition involves returning the quality of a product to a satisfactory state level (typically less than a virgin standard/new product) giving the resultant product a warranty less than of a newly manufactured equivalent.

[23] [START_REF] Paterson | End-of-Life decision tool with emphasis on Remanufacturing[END_REF] Repair is an activity of returning a used product in to "working order" by fixing/replacing specified faults in a product using service parts.

[22] [START_REF] King | Reducing waste: Repair, recondition, remanufacture or recycle?[END_REF][25] [START_REF] Krikke | Recovery Strategies and Reverse Logistic Network Design[END_REF] Recycle is an activity where discarded materials are collected, processed and used in the production of new materials or products.

[4] [START_REF] Ijomah | A Model-Based Definition of the Generic Remanufacturing Business Process[END_REF] ii

. Scoping of End-of-life Strategies

After defining potential EoL strategies, this step gives the decision maker an opportunity to take a look of defined product recovery strategies against a set of feasibility criteria for the refinement of viable EoL alternatives. The purpose of step 2 is to eliminate non-conforming scenarios during initial steps decision-making process based on various constraining influences such as technological, business, legislative and societal aspects that influence the feasibility of a particular EoL strategy. Table 2 shows a list of decision-making factors used in refinement of potential EoL strategies.

The screening process of the EoL strategies is mainly qualitative. The selected EoL options from the screening process will be considered in the following steps of the decisionmaking process. The selection of a potential EoL product recovery option should be based on the information available related to the activity and experience of the decision-maker [START_REF] Bufardi | Multicriteria decision-aid approach for product end-of-life alternative selection[END_REF].

Table 2: Categorization of EoL decision making factors

Category

List of key factors Ecological (Environmental) points (mPt) is used which is one-thousandth of a Pt. The end-of-life impact of a material for a specific strategy can be refereed from eco-indicator table [START_REF] Consultants | Eco-indicator 99 Manual for Designers[END_REF]. A positive point implies impact imposed on the environment while a negative impact infers impact which is avoided [START_REF] Lee | A multi-objective methodology for evaluating product end-of-life options and disassembly[END_REF].

Social indicator

Name: Number of employees to perform the scenario Unit: Integer number Goal: Maximizing number of employees

v. Analysis and Evaluation of End-of-life Strategies

This step involves the ranking of EoL strategies based on the information retrieved from step 2 and the selected environmental, economic and social indicators in step 3. The information and data gathered from each step is critically evaluated to select the most appropriate EoL treatment strategy [START_REF] Lamvik | The AEOLOS methodology[END_REF]. Due to the wide range of different multicriteria decision-making approaches, the choice of an appropriate method should be given great attention. It is critical for the decision maker to understand the problem, the feasible alternatives, conflicts between the criteria and level of uncertainty of the data before carrying out the choice to every multicriteria decision-making situation [START_REF] Bufardi | Multicriteria decision-aid approach for product end-of-life alternative selection[END_REF].

vi. Refinement of End-of-life Strategies and Final Evaluation

Once the analysis and ranking of potential EoL strategies is completed, further detail analysis should be applied by the decision maker to understand the consequences of selecting the best suitable strategy as a final solution. A critical evaluation of the potential best feasible product recovery strategy should be done against a set of criteria presented in table 2. This step may result in acknowledgement of the candidate strategy as a final solution or may lead to a new iteration of the approach. In case, the user found the result to be unsatisfactory, then the next EoL option is considered and evaluated in the same way as the previous candidate. Alternatively, the procedure will be repeated by considering a new set of EoL strategies and/or a new family of indicators [START_REF] Bufardi | Multicriteria decision-aid approach for product end-of-life alternative selection[END_REF] [START_REF] Lamvik | The AEOLOS methodology[END_REF].

End-of-Life Decision making factors

Findings from literature show that economic and environmental decision making factors are widely used to assess the viability of conducting remanufacturing while neglecting other equally important factors such as legislation and societal factors [START_REF] Goodall | A review of the state of the art in tools and techniques used to evaluate remanufacturing feasibility[END_REF] [START_REF] Doyle | Identifying the End of Life Decision Making Factors[END_REF]. Social decisionmaking factors are most valuable to provide feasibility analysis of adopting of remanufacturing at strategic level. Furthermore, there is lack of a holistic approach for analysing and evaluating different types of factors simultaneously.

Based upon a comprehensive literature review and feedback from expertise in the subject domain, key end-of-life decision-making factors used to assess the feasibility of remanufacturing for product recovery were identified and presented (see Table 2). The decision-making factors are categorized into business, technical, economic, environmental, legal and societal aspects.

Case study

To exemplify the application of product recovery multi-criteria decision tool (PR-MCDT), an illustrative example of an automotive engine is carried out to show how the approach can be used. At the end-of-life, an engine can follow different routes that have its own consequences from economic, environmental, societal and business point of view. In this specific case, a light fiat engine, is considered with the evaluation of its main components (cylinder block, cylinder head, pistons, connecting rods, crankshaft, Flywheel, Camshaft & Turbo) to simplify the complexity of the problem.

Step 1: Selection of Potential End-of-life strategies

The first step in this approach is to define the constitution of a set of potential EoL product recovery strategies. The selection of potential end-of-life strategies depends on the type of product and the associated product recovery option. In this specific case study, three potential end-of-life product recovery strategies are defined.

List of Potential Product Recovery Strategies

Alternative 1: Reusing the product with minor service cleaning, polishing)

Alternative 2: Remanufacturing Alternative 3: Recycling

Step 2: Scoping of EoL strategies In this step, potential EoL strategies are evaluated against list of criteria categorized in to legislative, technical, business and societal aspects which is presented in Table 2. Nonconforming scenarios will be eliminated from the list while the remaining ones will be evaluated in the following steps. The selection of relevant EoL strategies depends on the preferences of the user (recovery company), the objective of the problem, experience of the user and constraints from social, market, legislation and technology. It is assumed that potential EoL alternatives of the automotive engine fairly satisfies those requirements. In general, few EoL strategies are interesting for the decision maker from a list of potential recovery options.

Step 3: Selection of relevant indicator

Economic and environmental indicators are used for the evaluation of EoL alternatives. In this illustrative example, social indicator, number of employees to perform recycling and remanufacturing activities is assumed to be relatively matching. The evaluation of the EoL strategies with respect to the indicators is presented in Table 5. Before taking the final decision, EoL alternative 2 (Remanufacturing) should be examined in more detail with respect to the list of decision-making factors presented in table 4. Even if from a technical point of view, if remanufacturing of the automotive engine is possible, further investigation should be made to examine the selected strategy with list of pertinent decision-making factors like market demand and compliance with legislation. If it is realized that a the selected EoL option is unsatisfactory, another EoL option should be analysed again based on the ranking of the evaluation or the evaluation process is repeated with a consideration of alternative EoL strategies.

Conclusion

In this paper, we proposed a general product recovery multi-criteria decision tool (PR-MCDT) to evaluate the viability of remanufacturing at strategic level. The decision-making tool uses a holistic approach, under several often-conflicting criteria, to assess the feasibility of remanufacturing with respect to relevant business, legal, environmental, social and economic factors and by taking in to account the preferences of the decision maker.

Based on the analysis of literature and feedback form expertise, decision-making factors were also identified in regard to technical, economic, business, environmental and societal aspect. The paper also highlighted key decision-making criteria pertinent to consider in the decision-making process.

The paper dealt with important aspects related to the proposed approach such as definition of EoL strategies, selection of relevant indicators and exploitation of results. The proposed decision-making tool was also applied to an automotive engine case to illustrate the applicability of the approach. The results show that, remanufacturing is a feasible EoL option compared with repair and recycling strategies.

Figure 1 :

 1 Figure 1: Multicriteria Decision Tool (MCDT)

Table 1 : Product recovery strategies

 1

Table 5 : evaluation result of potential EoL strategies [32]

 5 Based on the evaluation of the potential EoL strategies with the relevant indicators, the decision-maker can eliminate potential options which have low result. In this illustrative example, it appears that the solution can be easily drawn based on table of evaluations. EoL alternative 2 (Remanufacturing) is the best compromise EoL strategy from an integrated point considering environmental, economic & societal indicators.

	EoL alternative 1	EoL alternative 2	EoL alternative 3
		Reuse	Remanufacture	Recycle
	Indicator 1(Economic)
	Revenue selling materials (€)			108.46
	Steel			5.74
	Cast iron			34.85
	Aluminium			68.44
	Revenue selling engine (€)	612.00	2562.00
	Operating costs	43.47	263.46	158.20
	Energy consumption	3.91	24.85	3.50
	Workforce costs	39.56	238.60	154.71
	Total Revenue (€)	568.53	2298.54	-49.75
	Indicator 2 (Env.impact)
	Kg CO2 eq. treatment process 561	10920	502
	Kg CO2 eq. recycling process			337154
	Kg CO2 eq. avoided	28978	72446
	remanufacturing
	Kg CO2 eq. avoided raw	640719	640719	640410
	material extraction
	Benefits	-669137	-702245	-302754

ICoR-17 -Linköping University (c) Authors of the paper

Acknowledgements

This research was fully funded by the "Marie-Sklodowska-Curie Innovative Training Network "Circ€uit"-, Circular European Economy Innovative Training Network", within Horizon 2020 Programme of the European Commission.

*Human health (HH) *Ecosystem Quality (EQ) *Resources (R) Legislation *Compliance with legislation/ EU legislation/WEEE *Compliance with new legislation Market *Customer demand (Market demand) *Competitive pressure

iii. Selection of relevant indicators

The implementation of EoL strategy to recover a product at its end-of-life has environmental, economic and societal impact. These impacts are measured by appropriate indicators to formulate a judgement on the selection of the best compromise EoL strategy. The selection of relevant indicators may be accomplished from a predefined list where the decision-maker decides based on EoL situation or develop his/her own individual indicators [16][19]. According to [START_REF] Bufardi | Multicriteria decision-aid approach for product end-of-life alternative selection[END_REF], indicators should be specified on the bases of direction of preference (maximization/minimization), scale (qualitative/quantitative) and unit of measurement. Table 3 shows an example of indicators.

iv. Assessment of End-of-life Strategies

Once the end-of-life indicators and potential product recovery strategies were selected, the next step will be an evaluation of each EoL strategies with respect to the defined indicators. Potential EoL alternatives (Alt 1, Alt 2, Alt 3 …Alt N) with respect to the evaluation criteria (I1, I2, I3 …IN) were presented in table 4 [START_REF] Bufardi | Multiple criteria decision aid for selecting the best product end of life scenario[END_REF]. After completing evaluation of strategies, strategies with a very bad (lowest) score will be eliminated. End-of-life options which do not fail to have on any indicator a worst value are considered on the second evaluation [START_REF] Bufardi | Multicriteria decision-aid approach for product end-of-life alternative selection[END_REF].

List of Indicators (I) Name


F Taia Alaoui

David Betaille

Valerie Renaudin

Pedestrian dead reckoning navigation with the help of A*-based routing graphs in large unconstrained spaces

Keywords: A* path-planning, Pedestrian Dead Reckoning, Map-aided, Particle filter

An A*-based routing graph is proposed to assist PDR indoor and outdoor navigation with handheld devices. Measurements are provided by inertial and magnetic sensors together with a GNSS receiver. The novelty of this work lies in providing a realistic motion support that mitigates the absence of obstacles and enables the calibration of the PDR model even in large spaces where GNSS signal is unavailable. This motion support is exploited for both predicting positions and updating them using a particle filter. The navigation network is used to correct for the gyro drift, to adjust the step length model and to assess heading misalignment between the pedestrian's walking direction and the pointing direction of the handheld device. Several data sets have been tested and results show that the proposed model ensures a seamless transition between outdoor and indoor environments and improves the positioning accuracy. The drift is almost cancelled thanks to heading correction in contrast with a drift of 8% for the non-aided PDR approach. The mean error of filtered positions ranges from 3 to 5m.

Introduction

Pedestrian Dead Reckoning (PDR) is widely adopted in the field of pedestrian navigation with handheld devices. It is particularly adapted to smartphone-based localization as inertial sensors can be designed in a MEMS (Micro-Electro Mechanical Sensors) technology, enabling them to be embedded in lightweight devices. Unlike GNSS receivers, inertial sensors are especially useful indoors as they allow standalone localization without sky visibility. Yet, due to gyro drift and step detection limitations, additional information is required to assist the PDR positioning process. For footmounted sensors, zero velocity update (ZUPT) calibration is exploited to adjust the positioning parameters by detecting stance phases within the gait cycle (static phase). Though, this calibration is not possible with handheld devices because of free hand motion and an increased difficulty to detect the stance phase. Outdoors, PDR can still be aided by GNSS [START_REF] Bancroft | Data Fusion Algorithms for Multiple Inertial Measurement Units[END_REF], but this is not feasible indoors because of signal unavailability and further measurements are needed. These could be provided by radio beacons or visual information. The first approach requires infrastructure deployment and training [START_REF] Hu | Design of Continuous Indoor Navigation System Based On INS and Wifi[END_REF], while the second necessitates a camera and further image processing for feature recognition [START_REF] Lee | An Indoor Positioning System Using Vision aided Advanced PDR technology without image DB and with motion recognition[END_REF]. A third possibility is to constrain the pedestrian's position using map information. Two main paradigms can be retained from previous work. Either walkable space is given by 2D maps delimited by obstacles [START_REF] Hafner | Evaluation of Smartphone-based Indoor Positioning Using Different Bayes Filters[END_REF], or it is given by a routing graph network that transforms the positioning process into a piecewise 1D model [START_REF] Willemsen | A topological approach with MEMS in smartphones based on routing-graph[END_REF]. In the first case, space is better explored but the map is not exploited further than for detecting static obstacles (e.g. walls). This means that no calibration is performed unless an obstacle is hit. On the contrary, routing graphs are much more constraining because the motion model is directly given by the graph network. Hence, their use has greater impact on the shape and accuracy of the trajectory and they have to be realistic enough to limit positioning errors. This paper focuses on routing graph-assisted PDR. In fact, routing graphs involve a simple motion model that allows both obstacle avoidance and the calibration of walking directions within straight line travels [START_REF] Willemsen | A topological approach with MEMS in smartphones based on routing-graph[END_REF]. Their use can even be extended to calibrating the step length model as reported in our previous work [START_REF] Alaoui | A multihypothesis particle filtering approach for pedestrian dead reckoning[END_REF]. Though, two major drawbacks make their use quite impractical and sometimes ineffective. First, there is no ubiquitous process that generates combined indooroutdoor pedestrian graph networks. Their construction can be time-consuming and sometimes inadequate with pedestrian motion. While the graph is expected to counterbalance PDR limitations, it may introduce additional errors due to inconsistencies with real displacement. Second, they do not handle "pseudo-random" trajectories within obstacle-free areas and during the transition between outdoor and indoor spaces as the freedom of motion increases. In these cases, routing graphs become inefficient as no motion assumptions can be relied on to design a realistic network. Grid-based models are well-suited for exploring large spaces [START_REF] Fetzer | On Monte Carlo smoothing in multi sensor indoor localisation[END_REF]. Yet, no calibration is possible with this approach. This even leads to overestimate the travel distance because linearity of movement is lost due to the grid structure. As a result, the PDR process is no longer assisted and the basic issues of gyro drift and uncalibrated travel distance not solved.

A solution to the vulnerability of routing graph assisted PDR navigation in obstacle free spaces is investigated in this study. The proposal is to make use of A* algorithm, which is commonly used for path planning within virtual worlds or for vehicle route guidance, in order to design a realistic routing graph in obstacle-free spaces. The graph is constructed on the basis of a set of waypoints that are crucial for pedestrian navigation. In fact, the latter are expected to balance the absence of geometrical constraints by providing strategic locations between which path computation with A* algorithm would be relevant. Therefore, the first contribution of this study is to mitigate the drift in the PDR approach even in large spaces. Moreover, the issue of accumulating errors during the transition between indoor and outdoor spaces is addressed by improving the routing graph relevance. This is accomplished through the computation of the likely paths a pedestrian may take from outdoor strategic locations to reach buildings entrance doors. In addition, the choice of A* algorithm is well-suited for handling different mobility profiles (e.g. personal disease that impacts the path of the target user) so that the approach can be customized for the specific needs of each pedestrian. Indeed, this is ensured thanks to a userrelative weighting of the map which is directly involved in the graph construction. This is because A* algorithm computes optimal paths according to walkability rates given by a map of navigation (e.g. walkability rate according to slopes, pavement, etc.). For instance, one would take the shortest path to cross an open area, whereas a disabled person may follow another path such as walking near walls or dedicated tracks.

A* pathfinding within navigation meshes

2.1. Overview of Navigation Mesh generation. A navigation mesh (NavMesh) is a set of several 2D or 3D polygons reachable by some given user [START_REF] Mendonça | Stealthy Path Planning Using Navigation Meshes[END_REF]. This structure can be obtained automatically using GIS software. Others make use of RECAST [START_REF]Recast: Detail API documentation for the members declared in Recast.h… • recastnavigation/recastnavigation@6f5c9f9[END_REF], an open source library within the community of virtual world designers. Several studies were conducted on how to construct navigation meshes from raw 3D models [START_REF] Afyouni | Spatial models for context-aware indoor navigation systems: A survey[END_REF], [START_REF] Van Toll | Navigation meshes for realistic multi-layered environments[END_REF]. NavMesh construction allows assigning weights that translate the polygon walkability rates according to some given characteristics (slope, type = {stairs, flat ground…}) and to the target user's mobility profile (ability to walk, to climb, to take the stairs, etc.). The weighting of the NavMesh is pretty important for considering the user's mobility profile in order to compute suitable routing paths. DETOUR [START_REF] Detour | GitHub[END_REF] is an example of open source software that allows A* pathfinding calculations on the basis of a NavMesh.

A shopping mall NavMesh generation.

The proposed method was applied to a shopping mall (Fig. 1). The latter is mapped in Google Maps with the names of stores, the mall principal gates and the set of walkable zones. We used QGIS 2.12.1 software to extract the map in the form of a raster with an 11cm resolution. Walkable zones were extracted by digitalizing the map and eliminating obstacles. Delaunay triangulation was used to create the NavMesh. All polygons were assigned the same weight as all experiments are conducted by healthy persons. This implies that A* computed paths will minimize distance and may result in a shape that can be obtained with different approaches [START_REF] Afyouni | Spatial models for context-aware indoor navigation systems: A survey[END_REF], [START_REF] Mortari | IGNM): A novel approach for deriving connectivity graphs for indoor navigation[END_REF].

A* pathfinding algorithm.

A* is an improved version of the Dijkstra algorithm [START_REF] Eraghi | Path length comparison in grid maps of planning algorithms: HCTNav, A #x2217; and Dijkstra[END_REF] . The latter aims at calculating the optimal path between two points according to one given cost function. Optimization applies to the travel distance, the time of travel, the expended energy, etc. This adaptability allows dealing with the issue of reduced mobility due to a disease or handicap. The A* optimization process is a discrete search scheme where space is modelled by a grid composed of cells. Each cell is explored according to an adjacency graph that models connections within the grid. Once visited, a cell is stored together with its assigned cost until the target cell is reached. The optimal path is given by the sequence of cells for which the sum of costs is minimal. A* is an extension to the principle of Dijkstra that calculates the cost of visited cells according to their distance from the starting point (which is the same as the Dijkstra cost) but also to their assumed distance to the target-cell in a heuristic approach. It is the a priori cost-to-go assumed lower than the actual (unknown) cost from the current cell to the target cell. A* is faster than Dijkstra algorithm due to the sorted exploration of cells according to their costs-to-go. The A* pathfinding scheme is applicable to any rasterized NavMesh where polygon weights are inherited by the raster cells.

Fig. 1: Walkable zones have been meshed (green polygons). Blank spaces are also walkable but the NavMesh is restricted to our zone of interest. The gray blocks are the stores considered as obstacles.

A* based routing graph generation

Routing graph generation on the basis of waypoints.

A waypoint is a punctual element that intervenes significantly in the pathfinding process. The significance of a waypoint is related to the role it plays in the pedestrian's decision making during her/his travel. Two situations are worth considering.

The first is when the pedestrian intends to reach a known destination within a familiar environment. In this case, it is obvious that (s)he follows the itinerary best suited to her/his mobility profile. Generally it is the shortest one but it could be any path depending on the cost-function that defines the user's mobility profile. This situation is handled by A* algorithm using a weighted NavMesh. The target destination is known if it is visible from the pedestrian's current position, for example if the pedestrian is located at a graph node that is visible from her/his target node. The second situation is when the pedestrian's destination is invisible due to obstacles. In this case, waypoints are built in order to discretize the possible target destinations which may be either final or intermediate destinations.

To better understand the implementation of waypoints, space has to be considered in relation to human spatial cognition. Indeed, walkable areas within a building floor could be corridors, rooms or halls. Outdoors, walkable areas are the space comprising sidewalks, footpaths, squares, car parks, etc. (Roads belong to the drivable areas.). This intuitive classification of space components allows defining scenes and extends what is called decision scenes to outdoor space and corridors. Decision scenes have been defined previously in [START_REF] Gaisbauer | Wayfinding Model For Pedestrian Navigation[END_REF] as the places that "can be entered and left and [are] physically bounded by buildings and other solid obstacles that prevent movement". The extension of decision scenes to corridors and footpaths in this paper is motivated by the fact that these elements have borders (physical borders for corridors and geometrical borders for footpaths) and that they can be accessed and left through specific points. For example two corridors are two separate scenes. They can be accessed and left through their intersection which is a building corner; a corridor and hall are also distinct scenes between which traffic flow is possible through their intersection. Two elements are important for pedestrian travel inside and between scenes. The first is related to the isovist, which is visible space from the pedestrian's current position. The isovist is mainly influenced by the presence of obstacles that have an impact on both visibility and the pedestrian's trajectory [START_REF] Victor | évaluation des déplacements piétons quotidiens[END_REF]. Second element is the set of portals that allow flows of pedestrians from scene to scene [START_REF] Gaisbauer | Wayfinding Model For Pedestrian Navigation[END_REF] such as corridor extremities, building gates, store entries or footpath extremities. Indeed, key elements in the pedestrian wayfinding behavior depend mainly on portals and obstacle borders (which determine the isovist). They materialize the fact that visibility and purpose (destination) are the parameters that give a shape to one's trajectory. In this study, the set of waypoints is composed of portals (corridor extremities, footpath extremities, outdoor/indoor doors) and obstacle corners. The graph construction based on waypoints is obviously more realistic than kernel-based methods where scenes are modelled by their centers [START_REF] Yang | Generation of navigation graphs for indoor space[END_REF], implying erroneously that the pedestrian walks systematically through the center of the scene. Besides, the structure of the graph is not entirely determined by the geometry of space but handles behavioral-based paths generated according to waypoints and reflecting pedestrian travel strategies.

Once waypoints have been constructed, they can be exported in a database that will assist the pathfinding process. Each pair of waypoints are related by a set of itineraries that are generated with A* algorithm and are part of the final routing graph. Multi-hypotheses motion is handled as in classical routing graphs by exploring several paths and keeping the ones that are best adapted to IMU measurements. The graph structure is stored in the form of a GIS database that contains the graph segments identifiers, their extremity nodes coordinates, their length and their connections within the graph in a node-connectivity approach. Each edge of the graph is oriented and has an entry node A and an exit node B. The node connectivity design is involved in the motion model used for filtering.

A* path planning for seamless transition between outdoor and indoor spaces and within obstacle-free areas.

Previous work demonstrated that indoor layout-based graphs are efficient for enhancing the PDR localization within geometrically constrained areas [START_REF] Willemsen | A topological approach with MEMS in smartphones based on routing-graph[END_REF], [START_REF] Alaoui | A multihypothesis particle filtering approach for pedestrian dead reckoning[END_REF]. Hence, the proposed method is applied only to the area that precedes the mall principal entrance as well as in the obstacle-free hall (Fig. 2). The latter represent the critical GNSS-deprived places where map information fails to provide a calibration to the PDR positioning process, emphasizing the main issue being addressed in this paper. Outdoors, the routing graph is constructed according to the geometry of sidewalks, crosswalks and footpaths.

Fig. 3 shows the generated routing graph within the test area. Three waypoints materialize the building entrance gate (drawn is black in Fig. 2). The focus is made on the waypoint where red paths intersect. It represents the left side of the mall entrance. Red paths are the A*-calculated itineraries relating the left side of the gate to other places of interest such as stores entries or the building's main corridors. The latter are modelled by a series of straight lines according to corridors main directions. Three paths relate the outdoor to the indoor space through the left side of the mall entrance.

This representation shows that straight line travels are privileged for practical and fast displacement. They are also more realistic in regards with pedestrian motion. In fact, some hypotheses are eliminated such as walking towards walls or making a series of turns to attend a place that can be reached straight ahead. Moreover, this provides a measurement of walking directions given by the graph segments in obstaclefree zones, which would be impossible if a grid-based approach had been applied, requiring additional information such as radio beacons signal to compensate for the PDR errors.

IMU fusion with GNSS and the routing graph with a particle filter

4.1.
Step detection and heading calculation. The step detection is realized after motion classification according to [START_REF] Susi | Motion Mode Recognition and Step Detection Algorithms for Mobile Phone Users[END_REF] by detecting peaks on the acceleration signal using an adaptive threshold algorithm [START_REF] V. Renaudin | Step Length Estimation Using Handheld Inertial Sensors[END_REF]. According to the same reference, a generic model is used to estimate the step length. This model relies on a set of three parameters trained on 10 subjects and is given by s = k. (h (a f + b) + c), (1) where s is the step length, h is the user's height, f is the step frequency, {a,b,c} the generic model parameters and k a scale factor that is expected to calibrate the model on the pedestrian.

Headings have been calculated with MAGYQ attitude estimation filter [START_REF] Renaudin | Magnetic, Acceleration Fields and Gyroscope Quaternion (MAGYQ)-Based Attitude Estimation with Smartphone Sensors for Indoor Pedestrian Navigation[END_REF] that fuses signals from a tri-axis accelerometer, a tri-axis gyroscope and a tri-axis magnetometer. Heading calculation considers different carrying modes such as the swinging mode or the texting mode. The device orientation in 3D-space is then obtained and the yaw angle deduced, giving the orientation of the device relative to the true North. 4.2. Calibration of the PDR parameters using GNSS positions and the routing graph. The PDR parameters to be calibrated are the step length and headings. The step length model needs to be adjusted to each user as the model parameters are dependent on the pedestrian's physiological features and gait cycle, whereas headings are potentially misaligned with the actual walking direction because of gyro drift and the device carrying mode. In fact, the device may be oriented towards a direction which is different from the walking direction. These errors are compensated by fusing the IMU with GNSS and the routing graph. The fusion is realized thanks to a particle filter that models the state (vector of unknown variables) by a set of particles (a set of sampled state vectors). The state vector contains necessary variables for determining the user's position and is presented in detail in section 4.3.

The step length model is adjusted to the pedestrian using both the graph and GNSS positions when available. In fact, the graph allows keeping the positions on plausible path(s). This directly impacts the travel distance. Besides, GNSS decreases the particles dispersion by bringing them next to the GNSS position. On the other hand, walking directions are given by path headings. They are particularly reliable indoors as the paths are calculated by the A*-algorithm or are directly given by the corridors main directions. The routing graphderived walking direction is then compared to the IMU pointing direction and the most likely path is selected. The difference between both headings gives the IMU angular misalignment with the pedestrian's walking direction.

The step length model is calibrated using a scale factor k . The latter is assumed constant as it depends mainly on the pedestrian's physiological parameters. It is modelled by a Gaussian variable with 1m mean and low variations as the step length is restricted to realistic values comprised between (0.60m -1m). k is expected to converge within straight line travels where an optimal calibration is guaranteed. Heading misalignment β is the angular difference between the actual walking direction and the pointing direction of the IMU. It is time-dependent because it varies with hand motion. Uniformly distributed samples of heading misalignment values are drawn with variations up to 30°.

Model implementation with a particle filter.

The state vector is given by { , , , , , }

p t Id Dp δ k β w ( 2
)
where Id is the edge identifier in the graph database, Dp is the curvilinear abscissa, δ represents the travel direction, which can be 1 if the edge is walked forward (from A to B) and -1 otherwise (from B to A), k is the scale factor that calibrates the step length model, β is the misalignment between the walking direction and the pointing direction of the handheld device, w is the particle weight, p and t refer respectively to the particle and time. Each particle is moved along the current path according to the travel distance given by the step length. The particle filter introduces white noise on each parameter to spread the particles and optimize the graph exploration. The particle position is determined by the identifier of the graph edge and its curvilinear abscissa. If the latter is above the edge length, one of its connected edges will be explored. The transition model is (3) where DU is compared to the current edge length L Id to compute the curvilinear abscissa Dp (Fig. 4), s is the step length, and G n is the modeling error represented by a Gaussian noise. α ∈ {0,1,2} and determines the number of steps to be made by one particle. It allows detecting over/mis-detected steps.

Fig. 4 gives the transition test process that allows computing the curvilinear abscissa Dp when an exit node B is exceeded (i.e. DU > 0). t is the time index. BA/ BB connectivity translates the search in the graph database in order to determine if the edge being explored is connected to the exceeded one via an A or a B node. The same logic is adopted if an entry node A is exceeded (i.e. DU <0).

Fig. 4: Dp computation according to DU

When GNSS positions are accurate enough, they are used together with the graph headings in order to calibrate both the step length model and PDR headings. The graph headings provide a measurement of walking directions which are compared to PDR headings, taking into account potential angular misalignment. This is performed according to

2 1 1 2 2 | ( ) | 1 . . ' 2 ( ) exp Id PDR PDR Id t t θ θ θ θ β σ σ w w - -     - +   = - + ∆ ∆       +     ∑ ( 4
)
P GPS P GPS E E N N -   ∆ =   -   (5) 1 1 / ( ² ²) 0 0 1 / ( ² ²) GPS GPS map E E map N N -   + Σ =   +     σ σ σ σ (6)
where t w is the particle weight, P E and GPS E are respectively the predicted East coordinate of a particle and its corresponding GPS observation. P N and GPS N are the same parameters along the North direction.

1 -∑ is the weighting matrix considering the graph and GPS position accuracies.

Indoors, some GNSS positions are still available but the latter are generally unreliable and are rejected. The rejection test is made by thresholding the signal to noise ratio (SNR) and the horizontal dilution of precision (HDOP). Only the graph headings are exploited in the weighting process. The weighting equation indoors is Once the particles have been weighted, some of them are assigned low weights and become useless in the process. Resampling is performed to duplicate the particles with high weights and delete the others. The particles amount is kept constant and their weights equal after each update in order to explore enough hypotheses of motion.

Experiments

Data collection with a HSGNSS and IMU in hand.

Three healthy volunteers (2 men, 1 woman) collected data with ULISS (Ubiquitous Localization unit with Inertial Sensors and Satellites) device (Fig. 5) held in hand (Fig. 6). Data were collected in both outdoor and indoor environments for two different device carrying modes. Two acquisitions were made in a texting mode and one in both swinging and texting modes. The average duration of acquisitions was about 14 minutes and the walking distance for each trial almost 1.5 km.

Different scenarios were chosen to perform the experiments. These will be explained in detail in section 6.1 where the text is enhanced by figures. [START_REF]ULISS[END_REF] comprises a 9 degrees of freedom inertial mobile unit, a high sensitivity GNSS (HSGNSS) receiver and antenna, a memory card and a battery. It delivers measurements that are time-stamped in GPS time. Inertial sensors and magnetometers provide measurements at a 200 Hz frequency.

Input data. ULISS device

The HSGNSS receiver operates in a standalone mode and delivers positions in real time at a 5 Hz frequency. Delivered positions are time-stamped in GPS time and have metric accuracies ranging from 2m up to 10m near buildings and tree shades. In this work, GNSS positions were interpolated at the step frequency in order to be fused with the PDR estimates of headings and step lengths. Besides collecting data with ULISS device in hand, all volunteers were equipped with an independent GNSS receiver carried in their backpacks and a small antenna attached to their caps. GNSS measurements were then post-processed in order to calculate reference trajectories by differential GNSS. This was performed using RTKLIB 2.4.2p12 software [START_REF]GitHub -tomojitakasu/RTKLIB_bin[END_REF]. Measurements from the embarked GNSS receiver and from a nearby base station were used to perform relative GNSS positioning. Obtained positions had decimetric accuracies up to several meters near buildings and other elevated features. Afterwards, thresholding was applied to three parameters in order to reject some outlier position estimates. These parameters are first the number of visible satellites, second, the ratio factor of ambiguity resolution and final, the horizontal dilution of precision. The resulting positions had precisions below 1m and were adequate for accuracy assessment in this work.

Results

6.1. Trajectory analysis. Fig. 7, 8 and 9 show the estimated trajectories (red) with the A*-generated routing graph and the particle filter described in section 4.3. The blue pattern corresponds to the PDR trajectory, while the green one is the reference trajectory obtained by differential GNSS. For both acquisitions, the travel distance seems to be overestimated. Yet, heading determination is more accurate for the first acquisition as the shape of the trajectory is more faithful to the building structure. The drift is higher for the second dataset and can be visually observed at the end of the PDR trajectory.

Where the drift is most important, the PDR trajectory presents major inconsistencies with the map. According the above figures, the drift has been corrected as the estimated trajectories are more compliant with the building structure and with footpaths in outdoor space. This has been achieved thanks to the proposed particle filter and to an increased conformity with pedestrian motion; demonstrating that the positioning accuracy for the texting scenario has been significantly enhanced using our approach. Unlike the two first acquisitions, there is a gap in the PDR trajectory because the subject did not go around the building as the two first subjects did. This gap is retrieved in the filtered trajectory.

Obviously, this scenario implies a less accurate heading determination and even an alteration in the travel distance estimation. In fact, the walking distance is underestimated for the swinging mode (First part of the travel until the subject reached the outdoor) and overestimated for the texting mode (Second part of the travel until the ending point). These errors can be noticed in the PDR trajectory.

The filtered trajectory shows that the drift has been corrected, resulting in a shape that is more compliant with the map and with the reference trajectory outdoors. Yet, the positioning accuracy seems to be decreased as compared with the two first datasets. Following section discusses the accuracy of estimated trajectories.

Error computation.

In order to assess the accuracy of our positioning method, filtered positions were compared to the reference positions interpolated at the step frequency. The average plane error ranges from 4 to 5 meters for the three datasets (Fig. 10). Accuracy is dependent on the quality of the PDR trajectory. Therefore, computed errors are more important for the second dataset considering the texting scenario. For the third acquisition that includes both swinging and texting, the accuracy is significantly decreased and more outliers (precisions above 15m) are detected due to mismatching errors (i.e. choosing wrong edges of the graph). These errors are mainly due to uncalibrated PDR parameters and are discussed in the following sections. For the first trial (Fig. 11), angular misalignment is comprised between -10° and +10° and varies around an approximate mean value of 0°. According to this distribution, the angular difference between walking directions and the pointing direction of the device is minimal. Hence, applied corrections compensate only for the gyro drift; which is rather logical regarding the texting mode scenario. For the second dataset (Fig. 12), estimated heading misalignment values are comprised between -15° and +12°. They are not equally distributed around 0° (e.g. between the 1 st and 2 nd minutes the mean value is over 5°). From this analysis, non-negligible hand motion can be assumed even if the subject intended to perform the experiment in the texting mode.

Fig. 12: Full texting scenario (2) Fig. 13 gives the estimated heading misalignment for the third trial. The first part (until 9.5min) of the travel was performed in the swinging mode. Hence, the estimated values vary significantly (-15° to +20°). Afterwards, oscillations between ± 5° over a mean value of 0° are retrieved, reflecting the texting mode of the acquisition. Step length model calibration. In this paper, a scale factor was introduced on the step length in order to calibrate the walking distance. Though, due to degraded GNSS signal and to short walking periods in outdoor space, the scale factor was not calibrated. In fact, regular and accurate GNSS positions are needed through straight line travels as cited in our introduction in order to calibrate the walking distance. These conditions were not verified during our experiments.

As a result, only corrected headings were relied on in the selection process. Therefore, distance calibration occurred only at junctions of the graph when a change in heading was detected. This explains the fact that the accuracy of filtered trajectories is still enhanced and compliance with the map improved.

Though, the uncalibrated walking distances caused some mismatching errors because direction change was detected too late or too early in the process. This can be noticed in the Northern part of Fig. 9. Indeed, while the pedestrian was intending to exit the building, the filtered trajectory indicates that he was walking towards a corridor. This happened for two reasons. First, the real trajectory is quite unusual in terms of pedestrian behavior. In fact, there is a change in heading (observed in the PDR trajectory) that is independent from space configuration, invalidating the assumptions that allowed constructing our navigation network. Second, the uncalibrated walking distance prevented the particles from reaching outdoor space at the right time. Another mismatching error occurred at the middle of the building (Between 300m and 400m North) because direction change was detected too late due to uncalibrated walking distance. Later on, the particle filter corrected for this error and converged over the right corridor thanks to the particle dispersion over the graph and to the multi-hypothesis approach.

Conclusions

A map-aided PDR approach where a routing graph is used as motion model has been proposed. Main contribution of this paper is A* algorithm adaptation to elaborate a pedestrian network that is capable of cancelling the gyro drift and the misalignment between the device orientation and the walking direction even in large spaces. These are GNSSdeprived and obstacle free areas where the limitations of mapaided PDR algorithms are most important. In fact, widespread map-aided PDR approaches do not compensate for these errors when pedestrian motion is unconstrained, mainly during the transition between outdoor and indoor spaces and when obstacles are absent. The A*-based routing graph mitigates the lack of obstacles thanks to a set of waypoints implemented according to human spatial cognition and to a weighted navigation mesh. This allows building a realistic motion model that meets the requirements of map-aided localization. Indeed, the proposed routing graph is well exploited because it gives prior knowledge about the pedestrian's destination and provides reliable measurements of walking directions. Results show that it is adequate for a seamless transition between outdoor and indoor environments and for enhancing the positioning accuracy even in large spaces. Achieved accuracies range from 3 to 5 meters and the drift is almost cancelled with the help of the routing graph. Though, some mismatching errors due to uncalibrated walking distance, especially while carrying the device in the swinging mode, might induce important positioning errors. Indeed, proper conditions of sky visibility and sufficient period of outdoor walking are prerequisite for the step length calibration before the pedestrian reaches indoor space.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Main area where A* algorithm has served for the routing graph construction

θ

  is the walking direction estimated by the PDR algorithms, Id θ is the predicted heading of the current path, β is the heading misalignment between the path direction and the PDR estimate, Id θ σ is the standard deviation of the path heading, and PDR θ σ is the standard deviation of the PDR walking direction estimate.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Ubiquitous Localization Unit with Inertial Sensors and Satellites

Fig. 7

 7 and 8 correspond to acquisitions performed in the texting mode. The starting point for both acquisitions lies on the top right extremity of the trajectories. Both subjects made a closed loop around the building, with an intermediate outdoor travel (Bottom side of the figures) before reaching the starting point back.

Fig. 7 :Fig. 8 :Fig. 9 :Fig. 9

 7899 Fig. 7: Estimated trajectory (in red). The Blue trajectory gives the PDR position estimates and the green one gives the reference trajectory (Only Texting)

Fig. 10 :

 10 Fig. 10: Plane errors for each dataset 6.3. Heading misalignment estimation. The figures below show the estimated heading misalignment values for each dataset.For the first trial (Fig.11), angular misalignment is comprised between -10° and +10° and varies around an approximate mean value of 0°. According to this distribution, the angular difference between walking directions and the pointing direction of the device is minimal. Hence, applied corrections compensate only for the gyro drift; which is rather logical regarding the texting mode scenario.

Fig. 11 :

 11 Fig.11: Full Texting scenario[START_REF] Bancroft | Data Fusion Algorithms for Multiple Inertial Measurement Units[END_REF]

Fig. 13 :

 13 Fig.13: Swinging-Texting scenario6.4. Step length model calibration. In this paper, a scale factor was introduced on the step length in order to calibrate the walking distance. Though, due to degraded GNSS signal and to short walking periods in outdoor space, the scale factor was not calibrated. In fact, regular and accurate GNSS positions are needed through straight line travels as cited in our introduction in order to calibrate the walking distance. These conditions were not verified during our experiments.As a result, only corrected headings were relied on in the selection process. Therefore, distance calibration occurred only at junctions of the graph when a change in heading was

Acknowledgements

This research work was done as part of the HappyHand (2015-2018) project which is funded by the 19th "Fonds Unique Interministériel" (French national R&D fundings).


Vincent Richard

Aurelio Cappozzo

Raphaël Dumas

Comparative assessment of knee joint models used in multi-body optimisation for soft tissue artefact compensation

Keywords: Multi-body optimisation, kinematic constraints, soft tissue artefact, degree of agreement, knee joint, human locomotion

Abstract

Estimating joint kinematics from skin-marker trajectories recorded using stereophotogrammetry is complicated by soft tissue artefact (STA), an inexorable source of error. One solution is to use a bone pose estimator based on multi-body optimisation (MBO) embedding joint constraints to compensate for STA. However, there is some debate over the effectiveness of this method. The present study aimed to quantitatively assess the degree of agreement between reference (i.e., artefact-free) knee joint kinematics and the same kinematics estimated using MBO embedding six different knee joint models. The following motor tasks were assessed: level walking, hopping, cutting, running, sit-to-stand, and step-up.

Reference knee kinematics was taken from pin-marker or biplane fluoroscopic data acquired concurrently with skin-marker data, made available by the respective authors. For each motor task, Bland-Altman analysis revealed that the performance of MBO varied according to the joint model used, with a wide discrepancy in results across degrees of freedom (DoFs), models and motor tasks (with a bias between -10.2° and 13.2° and between -10.2 mm and 7.2 mm, and with a confidence interval up to ±14.8° and ±11.1 mm, for rotation and displacement, respectively). It can be concluded that, while MBO might occasionally improve kinematics estimation, as implemented to date it does not represent a reliable solution to the STA issue.

Introduction

Joint kinematics estimation commonly relies on methods involving a mechanical model of the locomotor apparatus together with the stereophotogrammetric tracking of skin-marker trajectories. However, soft tissue artefact (STA), i.e., the relative movement between the skinmarkers and the underlying bones, introduces errors that jeopardise the information content of the skeletal motion estimation [START_REF] Leardini | Human movement analysis using stereophotogrammetry: Part 3. Soft tissue artifact assessment and compensation[END_REF][START_REF] Peters | Quantification of soft tissue artifact in lower limb human motion analysis: A systematic review[END_REF]. Since the artefact has a frequency content similar to that of bone movement, the problem cannot be solved by filtering [START_REF] Chiari | Human movement analysis using stereophotogrammetry: Part 2: Instrumental errors[END_REF].

Multi-body optimisation (MBO) is increasingly used with the intent to compensate for STAs.

The method embeds a rigid multi-body system and kinematic models of the joints involved, which means that the degrees of freedom (DoFs) of the joints are constrained [START_REF] Andersen | Kinematic analysis of overdeterminate biomechanical systems[END_REF][START_REF] Bonnechère | Physiologically corrected coupled motion during gait analysis using a model-based approach[END_REF][START_REF] Charlton | Repeatability of an optimised lower body model[END_REF][START_REF] Duprey | Influence of joint constraints on lower limb kinematics estimation from skin markers using global optimization[END_REF][START_REF] Lu | Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints[END_REF][START_REF] Ojeda | A method to evaluate human skeletal models using marker residuals and global optimization[END_REF][START_REF] Reinbolt | Determination of patient-specific multi-joint kinematic models through two-level optimization[END_REF].

Various mechanical linkages representing the knee joint and embedded in MBO have been described. These involve major simplifications with respect to real and subject-specific joints, and have less than six independent DoFs. The hinge joint [START_REF] Andersen | Kinematic analysis of overdeterminate biomechanical systems[END_REF][START_REF] Reinbolt | Determination of patient-specific multi-joint kinematic models through two-level optimization[END_REF] allows rotation about only the flexion-extension axis. The spherical joint, the most common representation of the knee in MBO [START_REF] Charlton | Repeatability of an optimised lower body model[END_REF][START_REF] Lu | Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints[END_REF][START_REF] Ojeda | A method to evaluate human skeletal models using marker residuals and global optimization[END_REF], allows all rotational movements but no translation. These models provide, in most cases, a rather inadequate 3D representation of the physiological movement of the knee [START_REF] Andersen | Do kinematic models reduce the effects of soft tissue artefacts in skin marker-based motion analysis? An in vivo study of knee kinematics[END_REF][START_REF] Clément | Can generic knee joint models improve the measurement of osteoarthritic knee kinematics during squatting activity?[END_REF]. Parallel mechanisms have also been used [START_REF] Duprey | Influence of joint constraints on lower limb kinematics estimation from skin markers using global optimization[END_REF][START_REF] Gasparutto | Validation of a multi-body optimization with knee kinematic models including ligament constraints[END_REF][START_REF] Valente | Effect of lower-limb joint models on subject-specific musculoskeletal models and simulations of daily motor activities[END_REF], the principle of which relies on compound joints representing an assembly of simple mechanical linkages. Although these models generally allow most rotations and translations, they couple the DoFs thereby prescribing displacements in a deterministic way (i.e., kinematics is imposed by the geometry of the model). A different approach consists in representing the behaviour of the knee directly by mathematically coupling the DoFs [START_REF] Bonnechère | Physiologically corrected coupled motion during gait analysis using a model-based approach[END_REF][START_REF] Li | The inaccuracy of surface-measured modelderived tibiofemoral kinematics[END_REF][START_REF] Scheys | Calculating gait kinematics using MR-based kinematic models[END_REF], with up to five DoFs driven by the flexion angle. A more recent modelling approach relies on a knee joint stiffness matrix and minimization of the relevant deformation energy [START_REF] Richard | Knee Kinematics Estimation Using Multi-Body Optimisation Embedding a Knee Joint Stiffness Matrix: A Feasibility Study[END_REF].

Based on a number of studies assessing MBO, it may be concluded that no fully satisfactory knee joint model has been found yet [START_REF] Andersen | Do kinematic models reduce the effects of soft tissue artefacts in skin marker-based motion analysis? An in vivo study of knee kinematics[END_REF][START_REF] Clément | Can generic knee joint models improve the measurement of osteoarthritic knee kinematics during squatting activity?[END_REF][START_REF] Gasparutto | Validation of a multi-body optimization with knee kinematic models including ligament constraints[END_REF][START_REF] Richard | Knee Kinematics Estimation Using Multi-Body Optimisation Embedding a Knee Joint Stiffness Matrix: A Feasibility Study[END_REF]. However, each of these assessment studies was performed on a single motor task (i.e., level walking, stepping-up, running or squatting). Moreover, some motor tasks (e.g., hopping, cutting) have not yet been investigated. Nor have all the abovementioned joint models been compared to date. Finally, existing comparisons have been based on the root mean square error between estimated and reference kinematics, without performing any deeper analysis regarding the relevant degree of agreement [START_REF] Mclaughlin | Testing agreement between a new method and the gold standard-How do we test?[END_REF].

This study aimed to comprehensively compare the performance of MBO embedding six different knee joint models selected from those proposed in the literature. This was made possible thanks to the availability of concurrently acquired reference, virtually artefact-free, bone kinematics and skin-marker data [START_REF] Cereatti | A proposal for the standardization of soft tissue artefact description and data exchange, and an open dataset on upper and lower limb artefacts during different motor tasks[END_REF]. The following motor tasks performed by able-bodied volunteers were analysed: level walking, hopping, cutting, running, sit-to-stand, and step-up. The degree of agreement between the reference and the MBOestimated joint kinematics was assessed by Bland-Altman analysis as well as using the relevant root mean square error and determination coefficient.

2.

Materials and Methods

MBO framework

In this study, each bony segment is fully located and oriented (i.e., bone pose) in the global reference coordinate system by means of natural coordinates [START_REF] De Jalon | Kinematic and Dynamic Simulation of Rigid and Flexible Systems with Fully Cartesian Coordinates[END_REF][START_REF] Dumas | 3D inverse dynamics in non-orthonormal segment coordinate system[END_REF]. Only the knee joint (i.e., the tibio-femoral joint) was considered in the study, meaning that the only segments involved in the MBO were the shank and the thigh.

Three types of constraints are typically used in MBO: driving constraints , rigid body constraints , and kinematic constraints . The constraints are split into two sets of equations. A set of "soft" constraints contains the equations that may be violated (i.e., ).

These constraints define the objective function of MBO, . A set of "hard" constraints contains the equations that must be fulfilled ( , ). In this framework, MBO is thus, to be regarded as a constrained optimisation problem. Note that a subset of the kinematic constraints , especially in the case of ligaments, may be considered as "soft" constraints and appended in the objective function, (with a weight matrix W).

The present study considered the following knee joint models described in the literature and implemented, using natural coordinates, as kinematic constraints within MBO:

-None: no joint model, where the relative movement of the tibia and the femur are independent from each other and joint dislocation is therefore possible (this is, of course, a borderline case of MBO);

-Spherical: spherical joint model, allowing the three rotations while impeding the three displacements;

-Hinge: hinge joint model, allowing only one rotation about the flexion-extension axis while impeding the other DoFs;

-Parallel: parallel mechanism with minimized ligament length variation, where two sphere-on-plane contacts stand for the contact between femoral condyles and tibial plateau.

m Φ r Φ k Φ m Φ   T mm f  ΦΦ r Φ k Φ 1 k Φ   11 T mm kk f              ΦΦ W ΦΦ
-Coupling: coupling curves between the DoFs, where internal rotation and adduction angles, as well as anterior and proximal displacements are functions of the extension angle through polynomial functions, and where lateral displacement is impeded; (1988). Note that the MBO embedding model None is actually equivalent to a single-body optimisation (e.g., [START_REF] Soderkvist | Determining the movements of the skeleton using wellconfigured markers[END_REF].

Joint kinematics estimation

Joint coordinate systems used to compute the kinematics of the knee joint were defined so as to satisfy the conventions for axes and Euler sequence proposed by the ISB [START_REF] Wu | ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion--part I: ankle, hip, and spine[END_REF].

The actual joint angles and displacements (extension, adduction, and internal rotation angles, and lateral, anterior, and proximal displacements) were computed from the natural coordinates [START_REF] Dumas | Joint and segment coordinate systems revisited[END_REF].

Experimental data

Right thigh and shank movement data from a single trial for each of the selected motor tasks, performed by able-bodied male subjects, were used for the analysis. These data were obtained from the datasets reported in [START_REF] Cereatti | A proposal for the standardization of soft tissue artefact description and data exchange, and an open dataset on upper and lower limb artefacts during different motor tasks[END_REF]. They included both virtually artefact-free bone-pose data, obtained using either pin-markers or biplane fluoroscopy, and concurrently acquired skin-marker data. The bony segment coordinate systems were defined based on bone k Φ anatomy and the reference positions of the skin-markers with respect to these coordinate systems were defined as their mean positions over the duration of the motor task. The data for level walking, hopping and cutting were from one volunteer (age: 22 years, mass: 63 kg, height: 1.75 m; [START_REF] Benoit | Effect of skin movement artifact on knee kinematics during gait and cutting motions measured in vivo[END_REF], while the data for running were from another volunteer whose anthropometric features were unknown [START_REF] Reinschmidt | Effect of skin movement on the analysis of skeletal knee joint motion during running[END_REF]. Relevant artefactfree data were obtained using pins inserted in the distal femur and proximal tibia. Data for the step-up and sit-to-stand tasks were from one male volunteer (age: unknown, mass: 83 kg, height: 1.75 m) and artefact-free data were obtained via biplane fluoroscopy [START_REF] Tsai | Effects of soft tissue artifacts on the calculated kinematics and kinetics of the knee during stair-ascent[END_REF]. Further details concerning the experimental set-ups, the definition of bony segment coordinate systems and relevant calibration and registration procedures used for the different datasets can be found in [START_REF] Cereatti | A proposal for the standardization of soft tissue artefact description and data exchange, and an open dataset on upper and lower limb artefacts during different motor tasks[END_REF] and in the above-mentioned references.

Assessment

For each of the six motor tasks, reference (i.e., artefact-free) femur and tibia pose and knee joint angles and displacements were reconstructed using pin-marker or biplane fluoroscopy data. Femur and tibia pose and knee joint kinematics were also estimated using the concurrently acquired skin-marker data and six MBO procedures each embedding one of the above-illustrated knee joint models (None, Spherical, Hinge, Parallel, Coupling and Elastic).

The degree of agreement between the joint angles and displacements derived from the six MBO procedures and the reference kinematics was assessed through Bland-Altman analysis [START_REF] Bland | Statistical Methods for Assessing Agreement Between Two Methods of Clinical Measurement[END_REF]. The bias (b) and confidence interval (CI; i.e., 1.96 standard deviation) were calculated. The root mean square error (RMSE) and coefficient of determination (R 2 ) were also calculated for the sake of comparison with previous studies.

Note that, when using the models Spherical and Hinge, displacements were null, thus the relevant coefficient of determination could not be computed.

Results

Full results are presented here for three motor tasks: level walking, hopping and cutting, while results for the other motor tasks are reported in Supplementary Material. Note that the results reported in the body of the paper and those in Supplementary Material lead to the same general conclusions.

Kinematics

Both reference joint angles and displacements and those estimated using MBO embedding the six joint models, are represented in Figure 1 FIGURES 1, 2 AND 3 ABOUT HERE

Statistics

The statistical parameters are reported in Table 1, Table 2 and Table 3 for level walking, cutting, and hopping, respectively. In most cases, the lowest bias b was obtained using MBO with model None, particularly for the joint angles. However, CI was higher than that obtained with the other models. As observed from the joint angle and displacement time histories, the different models embedded in the MBO provided disparate results. There were marked differences in the RMSE and R 2 values obtained from the different joint models, motor tasks and DoFs. TABLES 1, 2 AND 3 ABOUT HERE

Discussion

The aim of this study was to evaluate how effective MBO is in limiting the propagation of STA to joint kinematics estimates when various knee joint models are used. The joint models assessed in this study were selected from those proposed in the literature. Overall, the outcome was disappointing: no significant improvement in the quality of the estimates was provided by any of these joint models. For orthopaedic applications, estimates of knee joint kinematics should preferably be unbiased and it is desirable that they remain within a 2° and 2 mm confidence interval [START_REF] Iwaki | Tibiofemoral movement 1: the shapes and relative movements of the femur and tibia in the unloaded cadaver knee[END_REF]. The MBO method and the six joint models evaluated in the present study did not reach this level of accuracy.

Propagation of STA to joint kinematics

Given the fact that the most satisfactory outcomes were provided by the MBO method embedding model None, it can be inferred that most of the skin-marker datasets analysed entailed limited STA. Still, virtual joint dislocation, common with this joint model, is revealed by the high CI observed for joint displacements in most of the motor tasks. The results tend to support the conclusions presented in [START_REF] Andersen | Do kinematic models reduce the effects of soft tissue artefacts in skin marker-based motion analysis? An in vivo study of knee kinematics[END_REF], in that the use of MBO embedding mechanical joint models (Spherical and Hinge) did not eliminate or significantly reduce the effects of STA and did not improve the overall validity of knee joint kinematics derived from skin-markers placed on the thigh and shank. The conclusion is the same for those models that are regarded as more physiological (Parallel, Coupling and Elastic), although performance did depend on the motor task and the DoF considered. The latter models were developed as a better option with respect to Spherical and Hinge models because they allow the estimate of joint displacements. The good to excellent determination coefficients obtained for anterior and proximal displacements support the notion that they can provide physiological time histories. However, from the results on bias b, CI and RMSE, it is clear that joint displacement estimates remain inaccurate.

Most published evaluation studies could not count on the availability of artefact-free reference data and were based on skin-marker residuals [START_REF] Andersen | Kinematic analysis of overdeterminate biomechanical systems[END_REF][START_REF] Ojeda | A method to evaluate human skeletal models using marker residuals and global optimization[END_REF][START_REF] Reinbolt | Determination of patient-specific multi-joint kinematic models through two-level optimization[END_REF] or joint dislocation in the case of model None [START_REF] Lu | Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints[END_REF][START_REF] Ojeda | A method to evaluate human skeletal models using marker residuals and global optimization[END_REF][START_REF] Sholukha | Double-step registration of in vivo stereophotogrammetry with both in vitro 6-DOFs electrogoniometry and CT medical imaging[END_REF]. Only few studies provide a quantitative assessment of the effect of STA propagation to normal knee kinematics (i.e., non-pathological and nonoperated knees) using artefact-free reference data. For level walking and investigating six subjects with skin-and pin-mounted markers, [START_REF] Andersen | Do kinematic models reduce the effects of soft tissue artefacts in skin marker-based motion analysis? An in vivo study of knee kinematics[END_REF] reported a mean RMSE between 1.4° and 8.1° for extension, 1.5° and 9.7° for adduction, and 1.3° and 4.7° for internal rotation when using model Spherical. For ten subjects, performing a squatting task with reference markers mounted on an external fixator and using biplane radiography, Overall, the values obtained for RMSE (or b and CI) in the present study were lower than those reported in the literature for various motor tasks and with various means of assessing reference kinematics. This could be explained by the fact that the analysed datasets from skinmarkers involved small STAs, as demonstrated by the acceptable results obtained with model None. Assessing the degree of agreement through Bland-Altman [START_REF] Bland | Statistical Methods for Assessing Agreement Between Two Methods of Clinical Measurement[END_REF] allowed a deeper comparative analysis of the joint models. For instance, Hinge and Coupling models, with numerous "hard" constraints, typically resulted in a higher bias b for internal rotation. It was also expected that Parallel and Elastic models, which introduce "soft" constraints, would yield lower CI, but this did not occur systematically.

The contrasting results obtained using the same model for a given subject performing three different tasks show that a joint model may be better suited to one particular motor task, while other models might work better for other tasks. For instance, for hopping, model Coupling provided poor agreement on joint angles but good agreement on displacements, while model Elastic yielded opposite indications, with good agreement on angles and poor agreement on displacements. No doubt some assumptions on the dependent DoFs in the knee, as well as some physiological features (i.e., ligaments with minimal elongation, stiffness coefficients), all established from cadaveric measurements, are questionable when applied to certain motor tasks performed in vivo. In examining cutting and hopping, two motor tasks which have not been studied before, we selected a critical case to assess the joint models used in MBO. Knee kinematics during these tasks may actually be distinctive (i.e., large anterior and lateral displacements), due to muscle contractions as well as external loading. Moreover, the amplitude of STA may also be reduced due to the muscle contractions involved.

Limitations

This study presents a number of limitations. First, only one trial with one subject was studied for each motor task. Second, only two segments were considered in the MBO. Different results can be expected from using MBO that considers more segments and adding constraints for other joints (e.g., hip and ankle). As previously shown [START_REF] Duprey | Influence of joint constraints on lower limb kinematics estimation from skin markers using global optimization[END_REF], joint kinematics, when estimated using MBO, depends on the type and the number of joint models considered. Third, our study considered pin-marker data as a reference (i.e. artefact-free), and the drawbacks of this approach are well known [START_REF] Ramsey | Biomechanics of the knee: methodological considerations in the in vivo kinematic analysis of the tibiofemoral and patellofemoral joint[END_REF]. In particular, local anaesthesia of the soft tissues surrounding the zone where the pins were screwed into the bone may partially explain the small quantity of STA observed on the knee joint kinematics estimation provided by the MBO embedding model None. Fourth, modelling the lower limb with rigid segments and non-personalised joint constraints that impede or prescribe one or more DoFs may introduce errors. In that case the procedure, instead of compensating for STA as originally intended, adds modelling errors to it. Subject-specific joint models embedded in MBO, in particular models Spherical and Parallel [START_REF] Clément | Soft tissue artifact compensation in knee kinematics by multi-body optimization: Performance of subjectspecific knee joint models[END_REF], have been shown to have a beneficial impact on kinematics estimation.

Conclusion

The present study proposed a quantitative comparison of the estimation of knee joint kinematics using MBO embedding six different models versus reference (i.e., artefact-free) kinematics. As a conclusion, it came out that MBO cannot be considered a fully reliable solution to compensate for the STA. As a matter of fact, only invasive/ionising methods can currently give access to joint kinematics with the desirable level of accuracy (2° and 2 mm).

Nevertheless MBO remains a useful approach for further biomechanical computation such as inverse dynamics and musculoskeletal modelling. Coupling and Elastic considering the extension, adduction and internal rotation angles as well as the lateral, anterior and proximal displacements versus the reference. Except for R 2 (no unit), the units are degrees for angles and mm for displacements.

Figure and table captions

Extension Adduction Internal rotation b CI R 2 RMSE b CI R 2 RMSE b CI R 2 RMSE None
0,1 5,8 0,91 2,9 0,0 5,2 0,06 2,7 -0,2 4,1 0,21 2,1 Spherical -0,1 6,7 0,88 3,4 0,3 7,3 0,08 3,7 0,0 3,1 0,50 1,6 Hinge -2,0 6,7 0,93 3,9 3,2 2,3 0,01 3,4 7,9 3,7 0,00 8,1 Parallel -0,3 6,5 0,90 3,3 2,0 3,2 0,04 2,6 -0,9 3,7 0,60 2,1 Coupling 0,9 7,4 0,89 3,9 3,4 2,5 0,01 3,6 7,2 5,5 0,01 7,7 Elastic -0,2 6,2 0,90 3,2 1,7 4,5 0,11 2,9 -0,3 3,2 0,43 1,7

Lateral displacement Anterior displacement Proximal displacement b CI R 2 RMSE b CI R 2 RMSE b CI R 2 RMSE None
-0,2 7,2 0,02 3,7 0,9 10,1 0,05 5,2 0,1 14,8 0,05 7,5 Spherical 0,0 2,5 -1,3 0,0 3,0 -1,5 0,0 7,3 -3,7 Hinge 0,0 2,5 -1,3 0,0 3,0 -1,5 0,0 7,3 -3,7 Parallel -5,0 8,7 0,36 6,7 -1,6 4,4 0,15 2,7 -0,2 6,7 0,74 3,4 Coupling 0,0 2,5 0,01 1,3 -0,2 2,6 0,22 1,4 0,2 6,2 0,69 3,1 Elastic 0,6 6,3 0,31 3,3 -5,0 4,0 0,08 5,4 -1,2 9,7 0,13 5,0 0,0 1,3 0,98 0,7 -0,1 5,1 0,42 2,6 -0,1 4,8 0,26 2,5 Spherical 0,0 1,6 0,98 0,8 -0,2 3,6 0,40 1,9 -0,2 4,6 0,28 2,3 Hinge -1,1 2,5 0,96 1,7 -2,9 3,3 0,00 3,3 9,0 3,5 0,00 9,2 Parallel 0,4 1,9 0,97 1,1 -2,8 3,2 0,04 3,2 1,5 3,1 0,35 2,2 Coupling 1,7 3,6 0,89 2,5 -2,5 3,1 0,16 2,9 6,8 3,1 0,29 7,0 Elastic 0,0 1,6 0,98 0,8 0,1 3,1 0,44 1,6 -0,4 4,3 0,28 2,2

Lateral displacement Anterior displacement Proximal displacement b CI R 2 RMSE b CI R 2 RMSE b CI R 2 RMSE None
0,0 7,0 0,16 3,6 0,2 5,4 0,10 2,8 0,1 2,3 0,37 1,2 Spherical -0,1 1,6 -0,8 0,4 3,6 -1,9 0,0 2,5 -1,3 Hinge -0,1 1,6 -0,8 0,4 3,6 -1,9 0,0 2,5 -1,3 Parallel 1,4 1,7 0,33 1,6 -0,4 3,4 0,48 1,8 0,0 2,3 0,70 1,2 Coupling -0,1 1,6 0,00 0,8 -0,1 3,0 0,56 1,5 0,7 1,7 0,82 1,1 Elastic -0,7 4,1 0,26 2,2 -4,0 3,8 0,05 4,4 -2,6 2,5 0,05 2,9

-

  Elastic: elastic joint model based on the stiffness matrix, where all six DoFs are defined by the minimisation of the deformation energy. A detailed description of the MBO method embedding the different models (i.e., kinematic constraints ) can be found in Duprey et al. (2010), Gasparutto et al. (2015) and Richard et al. (2016). More specifically, for model Parallel, the model geometry was taken from Parenti-Castelli and Sancisi (2013) and, for model Coupling, the coupling curves between the DoFs were an adaptation, due to a different sign convention, of those provided by Walker et al.

  , Figure2, and Figure3for level walking, hopping, and cutting, respectively. Overall, the extension angles estimated through the six MBO procedures and the reference angles exhibited similar patterns. However, conflicting results emerged with regard to the other DoFs. Depending on the motor task and the DoF considered, large discrepancies were observed in the kinematic outcomes of the MBO for any given joint model. Note that the characteristics of the knee models meant that model Spherical provided null displacements, model Hinge provided null displacements as well as null adduction and internal rotation angles, and model Coupling provided null lateral displacement. No general tendency could be extracted from the overall analysis of the joint angle and displacement time histories. Nonetheless, the kinematics obtained using MBO embedding model None (i.e., no constraints) did not stand out in any way from the other five models.

For

  level walking, models Hinge and Coupling provided significantly higher b values for joint angles, with 7.9° and 7.2° respectively for internal rotation, while the lowest b value was provided by model None, with an absolute value between 0.0° and 0.2°. MBO embedding models Spherical, Hinge, Parallel, Coupling and Elastic generally provided lower CI than MBO embedding model None, except for extension. Models Parallel and Elastic provided the highest b values overall for displacements, with absolute values of up to 5 mm obtained for lateral displacement with model Parallel and for anterior displacement with model Elastic. For cutting, the lowest CI was provided by model Hinge for extension, by model Elastic for adduction, internal rotation and lateral displacement, and by model Coupling for anterior and proximal displacements, while for b the lowest values were obtained with model Elastic for extension, with model None for adduction, with model Spherical for internal rotation, and with models Spherical and Hinge for all displacements. For hopping, b was found null for extension with models None, Spherical and Elastic and for proximal displacement with Spherical, Hinge and Parallel. The lowest b values were obtained with models None and Elastic for adduction, with model None for internal rotation and lateral displacement, and with model Coupling for anterior displacement, while model None provided the highest CI values for adduction, internal rotation, lateral and anterior and the lowest for extension.

Figure 1 .

 1 Figure 1. Level walking. Knee joint angles and displacements estimated using MBO method

Figure 2 .

 2 Figure 2. Cutting. Knee joint angles and displacements estimated using MBO method

Figure 3 .

 3 Figure 3. Hopping. Knee joint angles and displacements estimated using MBO method

  None provided the highest RMSE, with an overall mean of about 2.3° and 3.6 mm for joint angles and displacements. Finally, for step-up with two subjects (skin markers and biplane fluoroscopy),[START_REF] Richard | Knee Kinematics Estimation Using Multi-Body Optimisation Embedding a Knee Joint Stiffness Matrix: A Feasibility Study[END_REF] reported RMSEs of up to 7.1° and 8.2 mm for joint angles and displacements with model None, up to 5.6° and 7 mm with model Spherical, up to 20.4° and 3.5 mm with model Parallel, and up to 6.4° and 3.6 mm with model Elastic. In the same study, the authors provide values for b and CI (maxima of the absolute values [b;CI]) of up to

	[1.9;13.9]° and [4.2;15.5] mm for joint angles and displacements with model None, up to
	[2.0;11.0]° and [5.3;9.2] mm with model Spherical, up to [19.4;12.7]° and [3.4;6.9] mm with
	model Parallel and up to [2.4;12.3]° and [2.4;7.0] mm with model Elastic.
	RMSEs for angles and displacements, respectively, were in the ranges [3.3° -8.1°] and [1.8
	mm -2.9 mm] (Spherical), [4.3° -9.2°] and [1.8 mm -2.9 mm] (Hinge), [4.5° -13.2°] and
	[1.8 mm -5.5 mm] (Coupling), and [3.2° -9.6°] and [1.6 mm -4.1 mm] (Parallel), at 64° of
	flexion. For running, with three subjects (skin-and pin-mounted markers), Gasparutto et al.
	(2015) found lower overall RMSEs with model Coupling for extension, adduction, and
	internal rotation (up to 2.0°, 3.2°, 5,5°, respectively) and for lateral, anterior and proximal
	displacements (up to 1.7 mm 4.5 mm and 2.5 mm, respectively). In the same study by

[START_REF] Clément | Can generic knee joint models improve the measurement of osteoarthritic knee kinematics during squatting activity?[END_REF]

found higher RMSEs, with values up to 21.8° for angles and up to 9.7 mm for displacements using knee models Spherical, Hinge, Coupling and Parallel. Mean

[START_REF] Gasparutto | Validation of a multi-body optimization with knee kinematic models including ligament constraints[END_REF]

, the mean RMSEs were below 2.5° and 4.1 mm and 2.2° and 3.1 mm with model Parallel and Spherical, respectively, for joint angles and displacements. Model

Table 1 .

 1 Level walking. Bland Altman bias (b) and confidence interval (CI) as well as correlation coefficient (R 2 ) and root mean square error (RMSE) when comparing results of the MBO method embedding the six different knee joint models None, Spherical, Hinge, Parallel, Coupling and Elastic considering the extension, adduction and internal rotation angles as well as the lateral, anterior and proximal displacements versus the reference. Except for R 2 (no unit), the units are degrees for angles and mm for displacements.

Table 2 .

 2 Cutting. Bland Altman bias (b) and confidence interval (CI) as well as correlation coefficient (R 2 ) and root mean square error (RMSE) when comparing results of the MBO method embedding the six different knee joint models None, Spherical, Hinge, Parallel, Coupling and Elastic considering the extension, adduction and internal rotation angles as well as the lateral, anterior and proximal displacements versus the reference. Except for R 2 (no unit), the units are degrees for angles and mm for displacements.

Table 3 .

 3 Hopping. Bland Altman bias (b) and confidence interval (CI) as well as correlation coefficient (R 2 ) and root mean square error (RMSE) when comparing results of the MBO method embedding the six different knee joint models None, Spherical, Hinge, Parallel,

Table 2

 2

		Extension		Adduction			Internal rotation
		b	CI R 2	RMSE	b	CI	R 2	RMSE	b	CI	R 2	RMSE
	None	-0,2	5,1 0,97 2,6	0,3	5,5	0,45 2,8	-0,6	4,4	0,41 2,3
	Spherical -0,3	4,8 0,97 2,4	0,4	6,0	0,33 3,1	-0,5	4,4	0,38 2,3
	Hinge	1,5	3,8 0,98 2,4	-10,2	6,6	0,02 10,7	7,3	5,6	0,01 7,8
	Parallel	3,6	4,9 0,98 4,3	-10,2	7,6	0,41 10,9	2,4	8,1	0,13 4,8
	Coupling 4,4	4,9 0,99 5,0	-8,4	5,1	0,56		1,5	11,1 0,32 5,8
	Elastic	0,1	4,5 0,97 2,3	-1,0	4,7	0,50 2,5	-0,6	4,0	0,52 2,1
		Lateral displacement		Anterior displacement		Proximal displacement
		b	CI R 2	RMSE	b	CI	R 2	RMSE	b	CI	R 2	RMSE
	None	-0,4	6,2 0,27 3,2	1,3	11,8 0,31 6,1	0,9	7,5	0,31 3,9
	Spherical 0,2	2,9 -	1,5	0,0	9,5	-	4,8	0,1	8,7	-	4,4
	Hinge	0,2	2,9 -	1,5	0,0	9,5	-	4,8	0,1	8,7	-	4,4
	Parallel	4,2	3,1 0,00 4,5	-2,4	8,6	0,78 5,0	1,0	7,7	0,93 4,0
	Coupling 0,2	2,9 0,01 1,5	-2,3	6,5	0,96 4,0	2,2	6,3	0,95 3,9
	Elastic	-2,1	2,3 0,40 2,4	-3,1	9,7	0,00 5,8	-4,5	8,6	0,03 6,3

Table 3

 3

	Extension		Adduction		Internal rotation
	b	CI R 2	RMSE	b	CI R 2	RMSE	b	CI R 2	RMSE
	None

Conflict of interest statement

No party having a direct interest in the results of the research supporting this article has or will confer a benefit on the author(s) or on any organization with which the author(s) is/are associated.