Gabriel Peyré Ceremade
email: gabriel.peyre@ceremade.dauphine.fr

Texture Synthesis and Modification with a Patch-valued Wavelet Transform

Keywords: Manifold of patches, non linear wavelet transform, image geometry, texture synthesis, texture modification

This paper models a texture as a 2D mapping onto a nonlinear manifold representing the local structures of the image. This manifold is learned from the set of local patches from an exemplar texture. A multiscale decomposition of this manifold valued representation is computed that mimics the orthogonal wavelet transform. The key ingredient of this decomposition is a geometric association field that drives the computations along the manifold. Iterated predictions leads to the computation of details coefficients over the features manifold. The resulting transform is invertible, non-linear and represents efficiently the local geometric structures of the exemplar. The multiscale coefficients of this transform are used to perform analysis and synthesis of textures.

Geometric Modeling of Images

Multiscale decompositions. The classical isotropic wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] cannot capture efficiently geometrical structures of textures. The local anisotropy of edges requires specific constructions such as the curvelets frame of Candès and Donoho [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise c2 singularities[END_REF] or the bandelets framework of Le Pennec and Mallat [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF].

Geometric structure propagation. The local geometry of images can be described as points on a curved manifold. Lee et al. [START_REF] Lee | The nonlinear statistics of high-contrast patches in natural images[END_REF] show how the edge manifold emerges from the set of patches extracted from natural images. To handle more complex texture features we propose to learn this manifold from an exemplar image. Computation are performed along this manifold using a local connexion describing how structures propagate in the image plane. Mumford first introduced this notion for edges propagation with the elastica model [START_REF] Mumford | Elastica and computer vision[END_REF]. Computer vision scientists such as Williams and Jacobs [START_REF] Williams | Stochastic completion fields: a neural model of illusory contour shape and salience[END_REF] proposed in some cases efficient approximations of this edge propagation field.

Texture synthesis. Pioneer work of Julesz [START_REF] Julesz | Visual pattern discrimination[END_REF] states filtering rules for the probabilistic characterization of textures. The wavelet domain modeling of Heeger and Bergen [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF] synthesizes cloudy textures.

Most successful approaches for texture synthesis in graphics are based on nonparametric copy of small patches from an original texture, see for example Efros and Leung [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF]. Recent approaches such as the method of Lefebvre and Hoppe in [START_REF] Lefebvre | Parallel controllable texture synthesis[END_REF] have re-casted this non-parametric patches sampling into a multiscale framework. A random shuffling of patches coordinates allows to deviate from the original texture. These modifications can be seen roughly as wavelet details that are added at various scales during the synthesis.

Our manifold description of textures paves the way between multiscale texture descriptions [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF] and more complex, non parametric sampling [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Lefebvre | Parallel controllable texture synthesis[END_REF].

Manifold valued analysis. Modeling data using a manifold structure can performed using non-linear estimation algorithms such as Isomap [START_REF] Tenenbaum | A Global Geometric Framework for Nonlinear Dimensionality Reduction[END_REF] or LLE [START_REF] Roweis | Nonlinear Dimensionality Reduction by Locally Linear Embedding[END_REF]. These dimensionality reduction procedures enable the analysis of libraries of images such as the ones studied by Donoho and Grimes [START_REF] Donoho | Image manifolds isometric to euclidean space[END_REF].

Processing of manifold valued functions is studied by Ur-Raman et al. [START_REF] Ur-Raman | Multiscale representations of manifold-valued data[END_REF] who propose a wavelet-like decomposition for such data. Our multiscale feature valued analysis presented in section 4 is inspired by this work but is adapted to geometric images and textures. In particular, our scheme takes into account the relative position of two features in the image, which enables the generation of anisotropic geometries such as elongated features.

Recent approaches to image synthesis in computed graphics use manifold modeling of textures. Matusik et al. define a manifold from a set of textures [START_REF] Matusik | Texture design using a simplicial complex of morphable textures[END_REF]. Lefebvre and Hoppe introduce a mapping of an image into a higher dimensional appearance space [START_REF] Lefebvre | Appearance-space texture synthesis[END_REF]. This embedding allows a synthesis with high fidelity and spacial variations.

Manifolds as Image Models

We model the local structure of textures using a set of manifolds {M w } w computed from an exemplar texture f e . The manifold M w is composed of patches of fixed size w that are extracted from the original texture f e

M w def. = p w t \ t ∈ [0, 1] 2 ⊂ L 2 ([-w/2, w/2] 2 ) ( 1
)
where

∀ x ∈ [-w/2, w/2] 2 , p w t (x) = f e (x -t). (2)
As explained by Grimes and Donoho [START_REF] Donoho | Image manifolds isometric to euclidean space[END_REF], the set M w might have a complex structure and be non-differentiable even for a simple images f e . Figure 1 shows some examples of patches extracted from a discretized texture. In numerical applications, the image f e is composed of n × n pixels, so M w is a set of patches considered as vectors of R m×m with m = wn.

Models and estimation.

A gray scale image is a mapping f : {0, . . . , n -1} 2 → R from pixels to real numbers. Our manifold description of images corresponds to a factorization

f : {0, . . . , n -1} 2 f -→ M w0 Φw 0 -→ R (3)
where the manifold M w0 is estimated from an exemplar f e (that might be different from f ) and the scale parameter w 0 > 0 represents the smallest size of a typical feature. Note that the mapping Φ w0 depends only on M w0 and thus on f e and can be computed independently of f . Factorization (3) extends to color images by considering each channel independently.

In this paper, we are interested by the mapping f which locates an image feature f (x) ∈ M w0 of size w 0 at location x in the image plane. The function Φ w0 controls the display of each feature, and we assume a simple sampling process

Φ w0 : p ∈ M w0 → p(0) ∈ R.
We note that recent works in harmonic analysis has focussed on the processing of such function Φ w0 ∈ L 2 (M) defined over a manifold domain M. In particular Fourier analysis of this function can be carried using eigenvector of the Laplace Beltrami operator. Szlam et al. [START_REF] Szlam | A general framework for adaptive regularization based on diffusion processes on graphs[END_REF] use this spectral decomposition to perform a non-linear filtering.

In this paper, instead of modifying Φ w0 , we propose to modify f which involves two key ingredients.

Estimating f such that factorization (3) hold at least approximately. To that end we perform a simple local best fit

f [k] def.
= argmin

g∈M w 0 nw0/2 t1=-nw0/2 nw0/2 t2=-nw0/2 |f [k + t] -g[t]| 2 . ( 4
)
The minimization of equation ( 4) extends to color images. In practice a color equalization is performed on f e prior to the matching as shown in image 2.

Performing signal processing for manifold-valued functions. The general frameworks of Ur Uaman et al. [START_REF] Ur-Raman | Multiscale representations of manifold-valued data[END_REF] is not suited for processing patches that have some spacial connexions. We develop our own computation model and derive a multiscale decomposition in this paper.

Examples of texture representations. In the case of a texture manifold M w0 learned from an exemplar f e using equation ( 1), the estimation of f from a given texture f defines a 2D-valued mapping

ϕ : [0, 1] 2 → [0, 1] 2 , where ∀ x ∈ [0, 1] 2 , f (x) = p w0 ϕ(x) ∈ M w0 ,
where p w t is the patch extracted from f e at some position t, as defined by equation (1). A similar coordinate mapping ϕ is generated by the texture synthesis algorithm of Lefebvre and Hoppe [START_REF] Lefebvre | Parallel controllable texture synthesis[END_REF][START_REF] Lefebvre | Appearance-space texture synthesis[END_REF].

Association Field over the Patch Manifold

The mean computation of two features g, g ′ ∈ M w must take into account the relative positions of the features in the image to process. The two features are thus assumed to be samples of some manifold valued function f such that f (x) = g ∈ M w and f (x ′ ) = g ′ ∈ M w . For the purpose of our multiscale processing framework, one only needs to compute such average for horizontal ω = (1, 0) and vertical ω = (0, 1) alignements of features. We thus assume that x -x ′ = wω. We describe the computation of the mean mapping in the horizontal case, the vertical one being similar. Starting from two feature elements g, g ′ ∈ M w , a new patch is computed by mixing the right side of g with the left side of

g ′ µ ω (g, g ′ )(x) def. = g(x 1 + w/2, x 2 ) if x 1 < 0, g ′ (x 1 -w/2, x 2 ) if x 1 0.
The mean m ω (g, g ′ ) ∈ M w of g and g ′ in the direction ω is then computed as the feature that best matches this mixing image ω (g, g ′ ) def. This numerical algorithm is a simple way to estimate an association field between feature elements. Similar ideas are used in graphics for texture synthesis using copy of patches, see for instance [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Lefebvre | Parallel controllable texture synthesis[END_REF]. More complex approaches relying on a biological or variational modeling of the association field could be used [START_REF] Mumford | Elastica and computer vision[END_REF][START_REF] Williams | Stochastic completion fields: a neural model of illusory contour shape and salience[END_REF].

= argmin

h∈M w ||µ ω (g, g ′ ) -h||.

Multiscale Feature-valued Transform

The multiscale transform decomposes a given patch valued image f : [0, 1] 2 → M w0 of n × n values into a set of n 2 coefficients {d j } J j=1 ∪ { fJ } that encodes at each scale j J def. = = log 2 (n 2 ) the details needed to reconstruct the features of the original texture. This transform is inspired by the interpolating manifold wavelet decomposition of Ur-Raman et al. [START_REF] Ur-Raman | Multiscale representations of manifold-valued data[END_REF] but with important differences that make it suitable for the analysis of local image features.

Manifold parameterization. In order to carry computations along the manifold M w we assume that it is globally parameterized by a bijective mapping

γ w : Ω ⊂ R p → M w . ( 5
)
Such a parameterization is hardly available in practice. In numerical applications, we compute an approximate mapping using the Locally Linear Embedding [START_REF] Roweis | Nonlinear Dimensionality Reduction by Locally Linear Embedding[END_REF].

The goal of this dimensionality reduction is to compute wavelet coefficients as in the parameter domain d = γ w (g) -γ w (g ′ ) ∈ R p of two patches g, g ′ ∈ M w . Ur-Raman et al. [START_REF] Ur-Raman | Multiscale representations of manifold-valued data[END_REF] define a wavelet difference d as a vector in the tangent plane of the that is tangent to the geodesic joining g and g ′ and whose magnitude is the length of the geodesic. We use a parametric approach to cope with our patch-valued data. The value of the dimensionality p should reflect the complexity of the local geometry of f e and we use p = 5 in our numerical applications.

Function [{d ω j }j,ω, fJ ] = multiscale transform fwd( f ).

Input : feature valued image f , Output : detail coefficients d ω j for j = 1, . . . , J and ω ∈ {(1, 0), (0, 1)}, coarse approximation at scale J, fJ .

Initialization : set f0 = f . For j = 1, . . . , J = log 2 (n), Set wj = 2 j /n.

(Expand) Consider fj[k] = fj-1[k] ∈ M w j .
For ω = (1, 0) and then ω = (0, 1) (Split) Compute even and odd set of coefficients: Reconstruction from the wavelet coefficients requires the computation of the inverse mapping γ -1 w . In numerical applications where M w is estimated using a discrete set of points, we use a nearest neighbor search

f e j [k1, k2] = fj[2k1, k2] and f o j [k1, k2] = fj[2k1 + 1, k2]. (Extract) Compute µ = m ω (f e j [k1, k2], f e j [k1 + 1, k2]) and set d ω j [k1, k2] = γw j (µ) -γw j (f o j [k1, k2
∀x ∈ R d , γ -1 w (x) def. = argmin g∈M ||γ w (g) -x|| 2 ℓ 2 .
Feature-valued Wavelet Decomposition In the following, we use cyclic boundary conditions, which means that we formally write v

[k] = v[k + n] for a vector v ∈ R n .
For the analysis of non-periodic textures, we use an extension by symmetrization along the boundaries of the image.

The algorithm starts with an initial feature valued mapping f0 = f of n × values. Each step generates successively coarser features valued images fj [k] ∈ M wj of n/2 j × n/2 j coefficients, with w j = 2 j /n.

Each column of the array fj is processed and we drop the column index to simplify the notations. The same computation is then performed on the rows of the array, as described in the pseudo code 1. The current vector fj is split into even and odd coefficients

∀ k = 0, . . . , n/2 j -1, f e j [k] def. = fj [2k] and f o j [k] def. = fj [2k + 1].
The even coefficients are used to predict the values of the odd coefficients using the averaging operator described in section 3

f o j [k] ≈ µ def. = m ω ( f e j [k], f e j [k + 1]) where ω = (1, 0).
The wavelet decomposition at scale j along direction ω encodes the prediction error, measured using differences computed over the parameter domain

d ω j [k] def. = γ wj (µ) -γ wj ( f o j [k]) ∈ R p .
The even coefficients are then copied fj = f e j in order to be further processed along the rows. Once fj as been processed, the resulting feature valued image of n/2 j+1 × n/2 j+1 values is copied into fj+1 considered as taking values into the coarser manifold M wj+1 . Note that at each iteration of the analysis, the assignment fj ← fj-1 assumes a switch from elements fj-1

[k] ∈ M wj-1 to elements fj [k] ∈ M wj .
The resulting wavelet transform can be inverted by retrieving the odd coefficients f o j from the available details d ω j and even coefficients f e j = fj

f o j [k] = γ -1 wj (γ wj ( f e j [k]) -d ω j [k]).
The pseudo code 1 implements the forward transform, whereas the pseudo code 2 implements the reverse transform. Figure 7 (b) shows the magnitude of the wavelet coefficients packed similarly to the orthogonal wavelet coefficients [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

Function f = multiscale transform bwd({d ω j }j,ω, fJ ). Input : detail coefficients d ω j for j = 1, . . . , J and ω ∈ {(1, 0), (0, 1)}, coarse approximation at scale J, fJ . Output : feature valued image f . Initialization : set fJ = f . For j = J = log (n), . . . , 1, Set wj = 2 j /n.

For ω = (0, 1) and then

ω = (1, 0) (Reconstruct) Compute µ = m ω (fj[k1, k2], fj[k1 + 1, k2]) and set f o j [k1, k2] = γ -1 w j (µ -d ω j [k1, k2]) and f e j [k1, k2] = fj[k1, k2]. (Join) Compute coefficients: fj[2k1, k2] = f e j [k1, k2] and fj[2k1 + 1, k2] = f o j [k1, k2]. (Swap) Switch lines and columns : fj[k1, k2] = fj[k2, k1]. End (Contract) Consider fj-1[k] = fj[k] ∈ M w j-1 . End Output f = f0.
Table 2: Pseudo code for the backward feature valued transform.

Texture Analysis and Synthesis

In this section, an input texture f e is used to provide an estimation of the feature manifolds M w for dyadic sizes w. The multiscale transform is used with these manifolds to perform texture synthesis and modification.

Texture Synthesis

An ideal interpolation of a coarse scale sets of coefficients can be constructed using the inverse multiscale transform. The user defines the four coefficients composing fj0 for the scale j 0 = 1 (they are chosen at random in our examples). The wavelet coefficients {d ω j } j,ω are set to zero in order to deviate as little as possible from the exemplar texture f e . Coefficients at scale j = 1 define four corners of the texture and successive steps in the backward transform (pseudocode 2) perform dyadic refinements which compute feature valued functions fj of increasing dyadic sizes. The corresponding images f j def.

= (Φ wj • fj ) have 2 j ×2 j pixels and can be seen as a non linear interpolation. The final synthesized texture is f def.

= (Φ w J • f ). The pixels of these intermediate images are sampled from the exemplar texture via a coordinate mapping f j [k] = f e [ϕ j [k]] and the mapping ϕ j is refined through the iterations. Figure 4 shows the iterations of the synthesis algorithm which refines f j and ϕ j .

Figure 5 shows two examples of texture synthesis. Note that the size of the synthesized textures can be arbitrary large. We use periodic boundary conditions for the backward transform so that the resulting texture is periodic (although the exemplar needs not be) and tiles the plane. In order to deviate from the original texture, one can use non-zero wavelet coefficients. More precisely, we define a parametric texture model over the wavelet coefficients using

∀ j > 0, ∀ ω, ∀ k ∈ {0, . . . , n/2 j -1} 2 d ω j [k] ∼ X(σ j ),
where ∼ means that each coefficient is drawn from a gaussian random variable X(σ j ) of mean 0 and variance σ j which allows a user control of the synthesis. Figure 6 shows different examples of texture synthesis using various spectral content σ j .

Texture Modification

In contrast to texture synthesis, texture modification requires both the forward and backward multiscale transform. The modification process takes as input an original texture f , which is to be modified according to some exemplar texture f e . This exemplar texture is used as the exemplar model for the feature valued multiscale transform.

In a pre-processing step, a feature valued mapping f is computed from f following equation (4). The forward transform described in listing 1 is used to decompose f as wavelets details {d ω j } j,ω . The texture can be modified in the wavelet domain by applying a non-linear thresholding that removes small Original Fig. 6. Texture generation with varying multiscale spectral distributions. The bottom line shows the spectral variation σj (coarse scales on the left). A small band indicate little deviation with respect to an ideal synthesis for the considered scale.

whose magnitude is bellow a threshold T

d ω j,T [n] def. = d ω j [n] if ||d ω j [n]|| T, 0 otherwise.
A more complex thresholding strategy could be used to include a normalization of the wavelet coefficients but we do not use it in our numerical experiments. Note that each wavelet detail d ω j [n] is a p-dimensional vector so its magnitude ||d ω j [n]|| is computed using the usual norm in R p . The thresholded wavelet coefficients d ω j,T are used to synthesize a modified feature valued function fT using the backward transform described in listing 2. The modified texture is computed as f T def.

= (Φ 2 J • fT ). On figure 7, bottom row, one can see how the wavelets coefficients are progressively shrinked toward zero. When T = max(d ω j ), the modification algorithm performs a pure texture synthesis since all wavelets coefficients are set to zero.

On figure 8 one can see how the modification process smoothly interpolates between the two textures by progressively adding and removing texture structures when wavelet coefficients are thresholded to zero.

Conclusion

In this paper, we have proposed a new texture model as a mapping into a non-linear features manifold. A multiscale decomposition of such a mapping is performed by using a pairwise association field. This non-linear analysis and synthesis framework allows to perform feature preserving texture generation and modification. This new model is promising for capturing multiple structures that exists in non-homogeneous textures using several target manifolds. Taking into account non-manifold constraints such as symmetries is another avenue for future work.

Fig. 1 .

 1 Fig. 1. Two examples of textures and some local geometric structures.

Fig. 2 .

 2 Fig.2. Example of coordinate mapping ϕ estimated for two different pairs of textures (fe, f ). The mapping ϕ(x) = (ϕ1(x), ϕ2(x)) ∈ R 2 is depicted using a red color for ϕ1 and green for ϕ2. Prior to computing ϕ, the colors of fe are equalized to match those of f . The reconstructed texture is Φw 0 • f .

Figure 2

 2 Figure2shows two examples of mappings ϕ estimated for target textures f of a size twice bigger (in pixels) than the size of the exemplar f e . For natural images, the mapping ϕ is usually piecewise linear and exhibit jump discontinuities. The right column of figure2shows the function Φ w0 • f ≈ f reconstructed from the manifold representation f . This reconstruction satisfies (Φ w0 • f )(x) = f e (ϕ(x)). If the learned manifold M w0 does not contains the local geometric structures of f , the reconstruction might differ from the original texture. On the example of figure 2 (upper row), one can see that the green line features of the pepper image are lacking.

Figure 3 Fig. 3 .

 33 Figure 3 illustrates the process of mean features computations M w .

Fig. 4 .Fig. 5 .

 45 Fig. 4. Progression of the synthesis algorithm.

Fig. 7 .

 7 Fig. 7. Progressive modification of a texture. (b,d,f,h) Magnitude S of the thresholded wavelets coefficients for T / max d ω j = 0, 0.1, 0.5 and 1 respectively. (a,c,e,g) Corresponding reconstructed textures.

Fig. 8 .

 8 Fig. 8. Two examples of texture modification. (a) Input texture f . (b) Exemplar texture fe. (c,d,e) Reconstruction from the thresholded wavelet coefficients for thresholds T / max(d ω j ) = 0.1, 0.5 and 1 respectively.

Table 1 :

 1 Pseudo code for the forward feature valued transform.

	])	and fj[k1, k2] = f e j [k1, k2].


Gabriel Peyré
email: gabriel.peyre@polytechnique.fr

Laurent D Cohen
email: cohen@ceremade.dauphine.fr

Landmark-based Geodesic Computation for Heuristically Driven Path Planning

No heuristic 1

Previous works

The extraction of shortest paths is a building block for a large class of applications ranging from graph theory to computer vision. The computation can be carried over a totally discrete setting (such as a graph) or over a discretized domain (such as an image or a 3D mesh).

Graphs and discrete computation. The canonical method to compute shortest paths on graphs is the Dijkstra algorithm, see for instance [START_REF] Cormen | Introduction to Algorithms[END_REF]. Fast exploration strategies have been used to speed up the computation. The A * algorithm [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF] makes use of an heuristic to reduce the search space. Other exploration strategies have been proposed in the field of artificial intelligence, such as such as IDA * [START_REF] Korf | Depth-first iterative-deepening: an optimal admissible tree search[END_REF].

Images and continuous setting. In order to extract geodesics for a continuous metric, the Fast Marching algorithm [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF] uses both a discrete scheme for derivation and a front propagation. A similar algorithm was also proposed in [START_REF] Tsitsiklis | Efficient Algorithms for Globally Optimal Trajectories[END_REF].

The minimal length properties of geodesics has been applied in computer vision, for example to solve global minimization problems for deformable models [START_REF] Cohen | Global Minimum for Active Contour models: A Minimal Path Approach[END_REF]. The continuous nature of this method is particularly attractive for image-based processings, for instance to extract tubular structures and centerlines in 3D medical data [START_REF] Deschamps | Fast Extraction of Minimal Paths in 3D Images and Applications to Virtual Endoscopy[END_REF].

Path planning: from discrete to continuous. Discrete computation on graphs give rise to numerous schemes to perform motion planification and the A * algorithm is intensively used for path-finding in video games [START_REF] Stout | Smart moves: Intelligent path-finding[END_REF]. For the case of an Euclidean metric, faster algorithms have been devised that exploit specific data structures such as visibility graphs [START_REF] Rohnert | Shortest paths in the plane with convex polygonal obstacles[END_REF].

The Fast Marching algorithm can be used to extract paths with a non Euclidean metric [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF]. The authors of [START_REF] Melchior | Consideration of obstacle danger level in path planning using A* and fast-marching optimisation: comparative study[END_REF] compare the Fast Marching and the A * algorithms to perform motion planification, but they do not propose a heuristic modification of the Fast Marching.

Shortest Path: Continuous and Discrete Algorithms

In this section we expose an unifying framework which includes the Fast Marching [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], Dijkstra [START_REF] Cormen | Introduction to Algorithms[END_REF] and A * [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF] algorithms together with our new heuristically driven propagation.

Front Propagation Methods for Shortest Path

We work on a discrete set of points and for each point x we have access to its neighbors y, defining the relation y ∼ x. These points can be embedded on a discretization grid (for the Fast Marching and for our method) or can be the vertices of a graph (for Dijkstra and A * ). Our goal is to compute the distance function U (x) = d(x 1 , x) to some starting point x 1 .

Initialization:

• Alive set: the starting point x 0 ; • Trial set: the neighbors of x 0 ; • Far: the set of all other grid points.

Loop:

• Let x be the Trial point with the smallest priority P(x); • Move it from the Trial to the Alive set;

• For each neighbor y of the current point x:

-if y is Far, then add it to Alive and compute a new value for U (y), -if y is Alive, recompute the value U (y), and update it if the new value is smaller, -recompute the priority P(y).

• If the end point x = x 1 is reached, stop the algorithm. The propagation methods label the points during the computation according to: • Alive is the set of points at which the distance value U has been computed and will not change; • Trial is the set of next grid points to be examined and for which an estimate of U has been computed; • Far is the set of all other grid points, for which there is not yet an estimate for U . Table 1 shows the main steps of the algorithms. Each algorithm must implement the following sub-functions

• A way to update the value U (y) at a given Trial point y.

This computation uses the values of U at the adjacent locations. This computation depends on the metric used, which can be defined on a graph (for discrete methods) or on the whole space (for continuous methods). We explain in sections 2.2 and 2.3 specific instantiations for the Dijkstra and Fast Marching algorithms. • A priority map P orders the set of Alive points according to some computational criterion. In the Fast Marching and Dijkstra algorithm, P(x) = U (x) is the current distance to the starting point. In our heuristical front propagation as in the A * algorithm, P(x) is chosen to minimize the number of visited points. We explain in section 4 how to actually construct a priority function P that makes use of a heuristic.

Discrete Case: Dijkstra

In the discrete setting, a symmetric weight g(x, y) is used to define the distance between two adjacent points x ∼ y of the graph. The length of a path

v = [v 1 , . . . , v m ], of ℓ(v) = m adjacent points v i is L(v) = g(v i , v i+1
), and we define the distance between two vertices

d(x 0 , x 1 ) = min m,v {L(v) \ ℓ(v) = m, v 1 = x 0 , v m = x 1 } .
The distance U (x) at a vertex x in the alive set is updated during the propagation according to

U (x) = min y∼x (U (y) + g(x, y)).
The shortest path v from x 0 to x 1 is tracked backward using v 0 = x 1 and v i+1 = argmin y∼vi U (y).

Continuous Case: Fast Marching

In R d , we are given a potential function g(x) > 0, and the weighted geodesic distance between two points x 0 , x 1 ∈ R d , is defined as

d(x 0 , x 1 ) = min γ 1 0 ||γ ′ (t)|| g(γ(t))dt , (1)
where γ is a piecewise regular curve with γ(0) = x 0 and γ(1) = x 1 . When g = 1, the integral in (1) corresponds to the length of the curve γ and d is the classical Euclidean distance.

The Fast Marching uses the fact that the function U satisfies the nonlinear Eikonal equation:

||∇U (x)|| = g(x). (2)
The distance U (x) = u at a point x = x i,j in the trial set is updated during the propagation according to the solution of

max(u -U (x i-1,j ), u -U (x i+1,j ), 0) 2 + max(u -U (x i,j-1 ), u -U (x i,j+1 ), 0) 2 = h 2 g(x i,j ) 2 .
This is a second order equation (the equation is written in R 2 for simplicity) and it can be solved as detailed for example in [START_REF] Cohen | Multiple Contour Finding and Perceptual Grouping Using Minimal Paths[END_REF].

The geodesic curve γ from x 0 to x 1 can be computed by extracting the parametric curve C(t) that solves the back propagation equation:

dC dt = - --→ ∇U with C(0) = x 1 .
This gradient descent is a local computation, and it only uses the value of U for a small fraction of the visited grid points. Note that these grid points are those located in the Alive set at the end of the front propagation procedure.

Heuristically Driven Front Propagation

In this section we explain our algorithm in the 2D setting, and show some numerical results that illustrate the main features of this method.

Propagation with a Heuristic

In order to minimize the number of Alive points at the end of the front propagation procedure, one should use a priority function P that tries to advance the front toward the goal point x 1 . In order to do so, we assume that, together with the the current weighted distance to the start point U (x) = d(x 0 , x), we have an estimate of the remaining weighted distance V (x) ≈ d(x 1 , x). Our heuristical front propagation algorithm follows the implementation of listing 1 with a priority map

P(x) = U (x) + V (x). (3)
The rationale behind the definition of P is that d(x 0 , x) + d(x 1 , x) is minimal and constant along the geodesic path joining x 0 and x 1 , see [9].

Algorithm A * . For discrete graphs, it has been shown that if the heuristic satisfies V (x) d(x 1 , x), then the extracted geodesic is a path of minimum length. This leads to the A * algorithm of [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF]. Various strategies have been proposed to devise admissible heuristics, see [START_REF] Goldberg | Computing the shortest path: A* search meets graph theory[END_REF] and the references therein.

Heuristically driven Fast Marching

On figure 2, one can see a front propagation, where we have used the oracle heuristic

V (x) = λ d(x 1 , x), with a parameter λ ∈ [0, 1].
The value λ = 0 corresponds to the classical Fast Marching propagation, which results in a large region of Alive points (colored in red). However, as we increase the value of λ toward 1, the explored region shrinks around the geodesic path that links x 0 to x 1 .

There are however two important issues with this ordering of the Trial set:

• This ordering can break the monotone condition that is required by the Fast Marching algorithm to produce a valid approximation of the continuous underlying distance function. We show in the numerical results presented in sub-section 4.2 that the geometric error on the extracted geodesic remains low. • We do not have an immediate access to the remaining distance d(x, x 1 ), since it would involve performing another front propagation from x 1 . We explain in section 4 how to overcome this problem.

Evaluation of an Heuristic

We cast the problem of finding a good heuristic into the problem of approximating the distance function d(x, y) between two points (x, y) by some function d. We define the heuristic using

V (x) = d(x 1 , x) ≈ d(x 1 , x).
This approximated distance d should satisfies d d in order not to deviate the propagation from the true shortest path (computed with V = 0). It must be fast to evaluate and can use a reasonable amount of pre-computed data.

To evaluate the quality of the heuristic, we propose two metrics • The approximation metric

E 1 (d, d) = |d(x, y) -d(x, y)| 2 dx dy,
is approximated during evaluation using a finite number of precomputed distance maps d(p i , x) to a set of points {p i } m i=1 chosen at random

E 1 (d, d) = 1 m m i=1 |d(p i , y) -d(p i , y)| 2 dy.
• The computation gain metric measures the usefulness of the heuristic for extracting geodesics between a set of couples of points {(x i 0 , x i 1 )} i . For each front propagation from x i 0 to x i 1 , we measure the area A i ( d) of the Alive set at the end of the propagation. We also measure the area A i (0) of the Alive set for a propagation without heuristic. We define the computation gain metric to be

E 2 (d, d) = 1 m m i=1 A i ( d)/A i (0).
When the geodesics queries are random, these two metrics tends to gives the same evaluation. However, the computation-gain metric is able to better adapt to typical applications where queries can be be highly non-uniform (for instance in road or tubular structure extraction).

No heuristic, 3D view of the potential

Landmark-based Front Propagation

In this section, we describe a particular implementation of the heuristic. This heuristic is computed using an approximated distance d evaluated with a set of pre-computed distances to landmark points.

Landmark-based heuristic

A new method for distance evaluation on graph has been introduced in [START_REF] Goldberg | Computing the shortest path: A* search meets graph theory[END_REF] as an admissible heuristic for the A * algorithm. We explain why this method can be useful for our heuristically driven Fast Marching and give a quantitative numerical study.

This method exploits the triangular inequality, since we have, for every couple of points x et y,

d(x, y) = sup z |d(x, z) -d(z, y)| .
This equality is still valid in the continuous framework of the Fast Marching, and in order to derive a useful heuristic, one can chose a small set z 1 , . . . , z n of Landmark points. The set of distance maps d k (x) = d(z k , x) is pre-computed once for all, and we define the approximation

d z1...zn (x, y) = sup k=1...n |d k (x) -d k (y)| . (4)
In the following, we drop the z 1 , . . . , z n dependencies and call the approximated distance d. We note that this approximation always satisfies the condition d d.

Figure 3 give an intuitive explanation of the efficiency of this approximation. In the ideal case, the geodesic γ k joining a landmark z k to x also pass through y. In this case, we have d(x, y) = d(x, y) and there is no approximation. But in most of the cases, this is not true, but a geodesic γ k passes close to y and d is indeed a good approximation of the real distance.

Landmark-based Path Planing

This landmark-based propagation can be used to speedup the computation of most path planing application, including path finding in video games or robot path planing.

γ γ k γ k x x z k y y y k Ideal case Real life case d ( x , y ) d ( y , z k ) - d ( x , z k ) = d ( x , y ) d ( y , z k ) - d ( x , z k ) > Figure 3
. Justification of the approximation properties of d.

On figure 4, one can see the active region explored by the propagation algorithm for an increasing number of landmark points. These points are chosen at random on the 2D image. The explored region progressively shrinks toward the central extracted curve when we add more points, since the heuristic is becoming more accurate.

On figure 5, one can see a numerical evaluation of the precision of the extracted geodesic. It reports the distortion ||γ -γ|| H between the original geodesic γ (computed without heuristic) and γ the one extracted with a heuristic. We use the symmetric mean-square Hausdorff metric

||γ -γ|| 2 H = 1 2 γ min y∈ γ ||x -y|| 2 2 dx + γ min x∈γ ||x -y|| 2 2 dx ,
since this captures well the geometric distortion caused by our partial propagation. On 2D maps containing important curvilinear features (such as the road on the example on the left), the distortion caused by the heuristic is nearly not noticeable. On the contrary, for relatively flat maps (such as the one depicted on the right) the distortion can be relatively high (of the order of a few pixels for a map of 256 × 256 pixels). This is because the salient features of this map catch the geodesic and avoid large lateral moves when the front propagation is modified.

In order to improve the quality of the heuristic and thus the resulting speed-up, one has to carefully seed the landmarks across the image. In the next section we set-up a complete framework for the evaluation of this sampling, together with strategies to get a high-quality seeding of the base points.

Seeding strategies for the landmarks

The quality of a sampling can be measured using the previously introduced metrics

E i (z 1 , . . . , z n ) = E i (d, d z1...zn ) for i = 1, 2.
In order to find sampling locations {z 1 , . . . , z n } for the landmarks, we use various strategies, among which : (a) A manual sampling that exploits specific knowledge of the 2D map g. This semi-automatic method is not studied in this paper. (b) A random sampling in the image. (c) A uniform sampling according to the distance d. This can be accomplished using the farthest point sampling procedure proposed in [START_REF] Peyré | Geodesic Remeshing Using Front Propagation[END_REF], where one chooses

z n+1 = argmax z ( argmin
1 k n d(z k , z)).
(d) A uniform seeding on the boundary of the domain.

(e) A greedy strategy minimizing the E 1 metric. We define recursively the locations

z n+1 = argmin z∈R E 1 (z 1 , . . . , z n , z),
where R is some small set of random candidates. (f) Same as (e) but using E 2 instead of E 1 .

Distance approximation evaluation.

A good seeding strategy should be able to put a landmark before any couple of points in the domain, near the geodesic path connecting these points. For a constant (i.e. Euclidean) metric, (d) is thus the optimal sampling scheme. However, for complex metric, this is no more the case.

On figure 6, one can see the approximated distance d(x 1 , x) to a target point x 1 , computed for n = 32 landmarks using various seeding strategies.

On figure 7, one can see a quantitative plot of the approximation error E 1 (x 1 , . . . , x n ). The error-driven seeding strategy (e) clearly performs the best as expected by definition. For a smooth potential g (left) the seeding strategy on the boundary (d) performs well, but it fails to capture the topology of complex 2D maps (right).

Computation-saving evaluation. In order to evaluate our algorithm in a real setting, we have set-up a complete framework using 500 typical queries {(x i 0 , x i 1 )} i for various 2D

Real distance to maps. Some resulting shortest paths queries are shown on the left of figure 8, they are relevant for many standard applications such as path finding in video games or robot path planing.

On the right of figure 8, we have shown the decreasing of the computation-saving metric E 2 (x 1 , . . . , x n ) for several sampling strategies. For the greedy sampling strategy (f), we have used another set of 100 typical queries, in order not to bias the computation of the error E 2 .

One can see that the computation-driven sampling strategy (f) clearly outperforms the other strategies. In particular, the approximation-based strategy (e) does not give good results, mainly because typical paths are not correlated to areas where the approximation error ||d -d|| is large.

Reducing Memory Usage

Cells Representation

Classical methods, such as using an octree data structure, can be used to reduce the memory usage of level set algorithms, for example in order to perform image segmentation [START_REF] Droske | An adaptive level set method for interactive segmentation of intracranial tumors[END_REF]. We choose to implement a simple data structure to reduce the memory usage by allocating the grid cell on the fly during the propagation. A typical cell data structure, for 2D computation, is: To be able to retrieve a given cell in constant time, we also store the list of allocated cells in a hash table. This is important because when a new cell is allocated, we need to connect it to the existing cells. This pointer-based representation of the neighboring relation is convenient to extract the geodesic with a gradient descent. There is some memory overhead due to the explicit storage of pointers to neighbors, but the fact that our scheme explores significantly less cells than the classical Fast Marching allows to save much more memory, as shown in next section. The computing time overhead due to the use of a hash table is about 40% in all our tests.

struct cell { // current

Computation from Compressed Data

To reduce the storage requirement of the distances to landmark d k (x), we implement a compression procedure. A typical implementation should • Allow random access of the value d k (x), without decompressing the whole data. • Be asymetric, since we do not care about the compression time, but we need fast access to distance values. • Give low decompression error, since we need an accurate heuristic and we need to satisfy as much as possible the condition d d. We use a two-steps procedure:

• Differential representation : A 1 (x) = d(x 1 , x) and A k+1 (x) = d(x k+1 , x) -d x1...x k (x k+1 , x),
where we use the approximated distance d x1...x k to give an estimate of d(x k+1 , x) using the available k first landmarks. This new representation remove the redundancy that exists between the different distance maps.

• We code each map A i using a vector quantization scheme [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]. In our tests, we use codeblocks of size 3 × 3, and we quantize the resulting vectors of R 9 using a codebook of size 256. This results in a memory gain of 36 : 1 with respect to storing single precision floats. The constraints we have on the computation from compressed data are similar to the process of rendering from compressed textures [START_REF] Beers | Rendering from compressed textures[END_REF]. We use the differential coding strategy to cope with the particular redundancy of distances functions. The compression produces neither noticeable artifact nor distortion on the extracted geodesics. The computation time overhead due to decompression is about 15% in CPU time.

Other Applications

Our landmark-based propagation can be applied verbatim to higher dimensional propagation and to triangulated surfaces in order to speed-up shortest path queries.

Constrained Path Planning

Geodesics can be used to compute the path of a robot with various shape and motion constraints [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF]. Basically, each additional degree of freedom add a new dimension to the domain in which the front propagation should be performed. Solving such high dimensional problems is time and memory consuming, so the use of a heuristic is highly desirable. In our experiment, the most important issue is the memory used by the full-grid classical Fast Marching, and the memory management strategy exposed in subsection 5 is crucial to scale to complex problems.

On figure 11, one can see two examples of path extractions in 2D with one rotational additional degree of freedom. This results in 3D front propagation, and the corresponding speed function is depicted on the left. Figure 9 shows the influence of the heuristic strength λ on the cells explored by the front propagation.

Geodesic Extraction on 3D Meshes

The Fast Marching algorithm has been extended to 3D meshes in [START_REF] Sethian | Computing Geodesic Paths on Manifolds[END_REF]. Our heuristically driven front propagation extends to 3D surfaces in a straightforward manner. We use a constant metric g = 1 and as the surfaces considered does

Conclusion

In this paper we have presented a simple modification of the Fast Marching to speed up the extraction of geodesics on images, higher dimensional data and triangulated surfaces. This modification takes into account a heuristics computed using a set of distances to landmarks. We reduce the memory requirement of the algorithm using specific data structures and a fast compression scheme.

Figure 2 .

 2 Figure 2. An example of 2D path planning. The set of alive points according to increasing heuristic is shown in red.

Figure 4 .Figure 5 .

 45 Figure 4. Example of propagation using a landmark-based heuristic with random seeding of base points.

(a) 1 x 1 x 1 xFigure 6 . 1 Figure 7 .

 111617 Figure 6. Graphical display of the approximated distance d(x1, x) for various seeding strategies.0

Typical shortest paths queries 10 logFigure 8 .

 108 Figure 8. Computation speed-up E2(z1, . . . , zn) obtained with a various number of landmarks. The areas colored in red correspond to the explored region without heuristic.

Figure 9 .

 9 Figure 9. Explored area for constrained path planning.

Figure 10 .

 10 Figure 10. Heuristically driven front propagation on 3D meshes, with an increasing number of landmark points.

Figure 11 .

 11 Figure 11. Examples of constrained path planning.not have boundary, the uniform seeding strategy (c) gives the best results for the E 1 metric. On figure12, one can see the approximated distance d(x 1 , x) for various number of landmarks. On figure1and 10, one can see the algorithm in action on various meshes, and for various number of landmark points.

Figure 12 .

 12 Figure 12. Approximation of the distance function d(x1, x) using d(x1, x) on a 3D mesh, with an increasing number of landmarks.

Table 1 :

 1 Pseudo-code for the common framework for front propagation.


Nicolas Reneaud

Emma Pierre

Raphael Zory

Frédéric Chorin

Luc Thomas

Pascale Chavet

Thelma Coyle

Eric Truchet

Stephane Puech

Matthieu Ollivier

Reneaud Nicolas
email: nicolas@tedorthopedics.com

Pierre Emma

Zory Raphaël

Chorin Frédéric

Thomas Luc

Chavet Pascale

Coyle Thelma

Truchet Eric

Puech Stephane

Ollivier Matthieu

Chabrand Patrick

Gerus Pauline

Affiliations

Validity of an instrumented knee brace compared to 3D motion navigation: A cadaveric investigation

Keywords: Connected knee brace, inertial measurement unit, joint kinematics, navigation, validity, reliability

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Background

Knee joint function can be altered by several pathologies, such as anterior cruciate ligament rupture (ACLr) and knee osteoarthritis (OA) [START_REF] García-Bógalo | Evaluation of intraoperative anterior cruciate ligament laxity reconstruction ଝ[END_REF][START_REF] Bytyqi | Gait knee kinematic alterations in medial osteoarthritis: Three dimensional assessment[END_REF][START_REF] Ren | Three dimensional knee kinematics and kinetics in ACL-deficient patients with and without medial meniscus posterior horn tear during level walking[END_REF][START_REF] Ismailidis | Kinematic changes in patients with severe knee osteoarthritis are a result of reduced walking speed rather than disease severity[END_REF]. A better understanding of knee pathologies has given rise recently to new knee diagnostic and analysis tools [START_REF] Paul | Five-repetition sit-to-stand[END_REF][START_REF] Ateef | Test-retest reliability and correlates of 6-minute walk test in patients with primary osteoarthritis of knees[END_REF][START_REF] Bernal | Reliability of the OptoGait portable photoelectric cell system for the quantification of spatial-temporal parameters of gait in young adults[END_REF]. Thus, motion capture and questionnaires have become the gold standard for practitioners for several pathologies [START_REF] García-Bógalo | Evaluation of intraoperative anterior cruciate ligament laxity reconstruction ଝ[END_REF][START_REF] Bytyqi | Gait knee kinematic alterations in medial osteoarthritis: Three dimensional assessment[END_REF][START_REF] Devita | Quadriceps-strengthening exercise and quadriceps and knee biomechanics during walking in knee osteoarthritis : A two-centre randomized controlled trial[END_REF]. The optoelectronic camera system (OCS) is the usual tool for quantifying body motion [START_REF] Bytyqi | Gait knee kinematic alterations in medial osteoarthritis: Three dimensional assessment[END_REF][START_REF] Gerbrands | Individual selection of gait retraining strategies is essential to optimally reduce medial knee load during gait[END_REF][START_REF] Favre | Gait analysis of patients with knee osteoarthritis highlights a pathological mechanical pathway and provides a basis for therapeutic interventions[END_REF][START_REF] Duffell | Gait adaptations with aging in healthy participants and people with knee-joint osteoarthritis[END_REF]. However, motion capture with OCS has few limitations. Firstly, the gait patterns can be affected by the laboratory environment compared to an ecological environment [START_REF] Huddleston | Ambulatory measurement of knee motion and physical activity: Preliminary evaluation of a smart activity monitor[END_REF][START_REF] De Müllenheim | Clinical Interest of Ambulatory Assessment of Physical Activity and Walking Capacity in Peripheral Artery Disease[END_REF][START_REF] Fusca | Validation of a wearable IMU system for gait analysis: Protocol and application to a new system[END_REF].

Secondly, the OCS is very expensive, restricted to a limited field of view in time and space [START_REF] Huddleston | Ambulatory measurement of knee motion and physical activity: Preliminary evaluation of a smart activity monitor[END_REF][START_REF] Fusca | Validation of a wearable IMU system for gait analysis: Protocol and application to a new system[END_REF] and time consuming. Finally, the OCS cannot be used in a home-based rehabilitation perspective.

To develop systems that overcome these disadvantages, recent research has focused on gait analysis using wearable sensors [START_REF] Tanimoto | Lower limb kinematics during the swing phase in patients with knee osteoarthritis measured using an inertial sensor[END_REF], and new systems have been developed to quantify human motion. For example, accelerometers are valid and reliable for the recording of activity, energy expenditure and acceleration data [START_REF] Poitras | Validity and reliability ofwearable sensors for joint angle estimation: A systematic review[END_REF]. Nevertheless, the precision of kinematic measurements is poor with an accelerometer sensor only due to the presence of drift during fast accelerations [START_REF] Dejnabadi | Estimation and visualization of sagittal kinematics of lower limbs orientation using body-fixed sensors[END_REF][START_REF] Antonello | Acceleration Measurement Drift Rejection in Motion Control Systems by Augmented-State Kinematic Kalman Filter[END_REF]. Combining a gyroscope with an accelerometer may solve this problem. Indeed, gyroscopes can estimate sensor orientation by integrating the angular velocity on the X, Y, Z axes relative to the sensor. Yet this sensor is subject to a drift in angular velocity due to signal integration, a well-known problem in robotics and motion analysis. By mixing the complementary data from accelerometers and gyroscopes, the two drifts can be minimized with a fusion algorithm like a Kalman filter [START_REF] Poitras | Validity and reliability ofwearable sensors for joint angle estimation: A systematic review[END_REF][START_REF] Šlajpah | Kinematics based sensory fusion for wearable motion assessment in human walking[END_REF][START_REF] Valenti | Keeping a good attitude: A quaternion-based orientation filter for IMUs and MARGs[END_REF]. The association of the two sensors is called an inertial measurement unit (IMU). Body sensor networks with IMUs such as Xsens MVN (Xsens Technologies BV, Enschede, NL) [START_REF] Nazarahari | Sensor-to-body calibration procedure for clinical motion analysis of lower limb using magnetic and inertial measurement units[END_REF] or Perception NEURON Pro (Noitom, Miami, FL, USA) [START_REF] Sers | Validity of the Perception Neuron inertial motion capture system for upper body motion analysis[END_REF] provide quantitative movement analysis in ecological environments.

In a home-based rehabilitation perspective, the knee motion analysis has become necessary to supervise patients through biofeedback. To deal with that, a new instrumented knee brace with two IMUs was developed. These sensors calculated the movement of the thigh relative to the shank. Moreover, this integrated system compares to simple IMU sensors described previously did not necessitate clinicians or physiotherapist for calibration or anatomical placement due to its ease of use. Indeed, once the connected knee brace is in place, it is possible to quantify knee motion in an ecological environment to reduce bias in laboratory conditions with a smartphone [START_REF] Lindemann | Spatiotemporal gait analysis of older persons in clinical practice and research: Which parameters are relevant?[END_REF][START_REF] Stellmann | Ecological validity of walking capacity tests in multiple sclerosis[END_REF][START_REF] Mcfadyen | Using dual task walking as an aid to assess executive dysfunction ecologically in neurological populations: A narrative review[END_REF]. Palpation to find anatomical landmarks is not required, nor is it necessary to fix markers or IMUs with different systems, as it is the case with other systems [START_REF] Bytyqi | Gait knee kinematic alterations in medial osteoarthritis: Three dimensional assessment[END_REF][START_REF] Clément | Soft tissue artifact compensation in knee kinematics by multi-body optimization: Performance of subject-specific knee joint models[END_REF]. Moreover, this technology offers free use of a workstation for motion computation through cloud computing on the smartphone app, and biofeedback can be displayed in real time with smartphone computation. All these data allow a quantified rehabilitation for physiotherapists or clinicians via success scores and pain scale charts. Then, this system enhance physical rehabilitation by motivating the user with visual data of progression and informative videos [START_REF] Leardini | Validation of the angular measurements of a new inertial-measurement-unit based rehabilitation system: Comparison with state-of-the-art gait analysis[END_REF].

However, this knee brace, like IMU systems or OCS, is sensitive to soft tissue artifacts (STA), a problem inherent to all skin-mounted tracking system. A systematic review noted the effects of STA in the lower limb, finding that several factors influence the results, including sensor orientation [START_REF] Peters | Quantification of soft tissue artifact in lower limb human motion analysis: A systematic review[END_REF]. STA was found to reach more than 30 mm of magnitude on the thigh segment and more than 15 mm on the tibia [START_REF] Peters | Quantification of soft tissue artifact in lower limb human motion analysis: A systematic review[END_REF]. To overcome this problem, OCS investigation with bony pins would be highly accurate [START_REF] Lafortune | Three-dimensional kinematics of the human knee during walking[END_REF][START_REF] Benoit | Surface marker cluster translation, rotation, scaling and deformation: Their contribution to soft tissue artefact and impact on knee joint kinematics[END_REF], but this procedure is unethical and should not be applied to living patients. Along the same lines, the fluoroscopic 3D-matching technique is invasive due to radiation exposure [START_REF] Tsai | Quantification of three-dimensional movement of skin markers relative to the underlying bones during functional activities[END_REF]. Recently, several studies have used commercial navigation devices to assess knee kinematics during surgical procedures for ACL reconstruction or total knee arthroplasty [START_REF] Imbert | Knee laxity modifications after ACL rupture and surgical intraand extraarticular reconstructions : intraoperative measures in reconstructed and healthy knees, Knee Surgery[END_REF][START_REF] Grassi | ACL reconstruction with lateral plasty reduces translational and rotatory laxity compared to anatomical single bundle and nonanatomical double bundle surgery: An in vivo kinematic evaluation with navigation system[END_REF][START_REF] Feczko | Tibial component rotation in total knee arthroplasty[END_REF]. Navigation is composed of three parts: computer platform, tracking system, and rigid body marker. Tracking system visualizes the rigid body markers through OCS and calculates motion in real time with the tracking systems. Markers are fixed to the patient's bones surgical instruments to track joint motion or target objects after a referencing procedure which consists in determining bone positions compared to marker positions [START_REF] Bae | Computer Assisted Navigation in Knee Arthroplasty[END_REF]. For example, this system via bony pins and OCS is a computer-assisted surgical navigation device to help tunnel placement during ACL reconstruction and can quantify knee laxity and validate the surgical procedure [START_REF] Grassi | ACL reconstruction with lateral plasty reduces translational and rotatory laxity compared to anatomical single bundle and nonanatomical double bundle surgery: An in vivo kinematic evaluation with navigation system[END_REF][START_REF] Nakamae | Evaluation of the intraoperative kinematics during double-bundle anterior cruciate ligament reconstruction using a navigation system[END_REF].

A more ethical solution might be considered to begin with cadaveric investigations. Indeed, tibiofemoral motion in the unloaded cadaveric knee is almost the same as that in the living knee [START_REF] Hill | Tibiofemoral movement 2: The loaded and unloaded living knee studied by MRI[END_REF][START_REF] Worlicek | The influence of varus and valgus deviation on patellar kinematics in healthy knees: An exploratory cadaver study[END_REF].

Given the innovative aspect of the knee brace, its accuracy and reliability must be confirmed before the system can be used more extensively. In the present study, we thus investigated the accuracy and reliability of this new system of knee motion quantification and compared it to the gold standard during specific movements well-known to physiotherapists and orthopedic surgeons.

Material and methods

Subjects

Thirteen intact knees from seven frozen cadaveric males were used for this study. The specimens were thawed at room temperature for at least three hours and they showed no sign of degeneration. Exclusion criteria were signs of knee instability (Lachman test) and signs of surgical procedures on knee or hip, evidence of a knee or hip prosthesis, ACL reconstruction, etc. The mean age of the specimens was approximated at 85 years-old approximately due to a lack of anthropological information.

Three clusters of four reflective markers were fixed by one screw on the pelvis and two screws on the femur and tibia to record the knee kinematics. The femoral and tibial screws were fixed in the bone diaphysis outside the joint capsule with a 2-cm incision to avoid influencing the soft tissue. After fixing the three clusters, seven anatomical landmarks were identified with reflective markers fixed onto the skin with cyanoacrylate glue. The medial malleolus (MM), the lateral malleolus (LM), the medial and lateral epicondyle of the femur (MF, LF), the medial and lateral epicondyle of the tibia (MT, LT) and the greater trochanter (GT) were located.

Kinematic data were simultaneously obtained with seven cameras (Vicon Motion Systems Ltd.

Oxford, UK): four MX-T20 and three MX-T40 and with the instrumented knee brace with a sampling frequency of 100 Hz. The instrumented knee brace is composed with two IMU ICM-20948. The knee sleeve was designed to reduce compression around the knee by using soft tissue in opposition to medical knee sleeve designed to change the knee joint kinematic [START_REF] Mohd Sharif | Biomechanical and functional efficacy of knee sleeves: A literature review[END_REF].

The knee sleeve compression around the knee was minimal to avoid change in knee joint movement but enough to limit the movement of the IMU sensors. IMU's knee brace had a 16 bits resolution for accelerometers and gyroscope. However due to Bluetooth Low Energy communication, these resolutions were downgrade to 10 bits with respectively, a sensitivity of 7.81 mg/LSB for a ±4 g range and a sensitivity of 0.98 dps/LSB for a range of ±500 dps.

Experimental protocol

To evaluate the position of the clusters relative to the bones, a functional movement protocol was generated with the leg extended on the table in a static position then with the leg fully extended realizing hip circular movement. The hip joint can be modeled as a ball-and-socket joint, and the hip joint center (HJC) location is characterized by a point invariant in any position of the joint [START_REF] Jenny | Reproducibility of intra-operative measurement of the mechanical axes of the lower limb during total knee replacement with a non-image-based navigation system[END_REF]. The coordinates of the femoral head location were obtained by an optimization method to minimize the HJC movement in the circle movements compared to the hip reference.

The distal femur location was determined at the middle of the MF and LF markers. The proximal and distal tibia locations were calculated between the MT and LT markers and the MM and LM markers, respectively.

The mechanical axis of the femur was defined by the line containing the femoral head location and the distal femur location. The femoral condylar axis was defined by the MF and LF markers.

The femoral frontal plane was defined by the plane containing the mechanical axis and the femoral condylar axis. The femoral sagittal plane was defined with the femoral axis and the cross-product of the femoral axis and the femoral condylar axis [START_REF] Jenny | Reproducibility of intra-operative measurement of the mechanical axes of the lower limb during total knee replacement with a non-image-based navigation system[END_REF]. The mechanical axis of the tibia was defined by the line containing the proximal and distal locations of the tibia. The tibial condylar axis was defined by the MT and LT markers. The tibial frontal plane was defined by the plane containing the mechanical axis and the tibial condylar axis. The tibial sagittal plane was defined with the tibial axis and the cross-product of the tibial axis and the tibial condylar axis. To estimate the orientation of the two IMUs of the instrumented knee brace, we used a procedure implemented in a smartphone application.

The operator performed four movements, five times each: a flexion with a full extension to 100° of flexion (FLEX100), a maximal internal/external rotation with a flexion of 30° (ROT30), a pivot shift test (PS) and a flexion from full extension to maximal knee flexion under constraint (FLEXMAX).

Calculation

Analysis

The synchronization of the two systems was performed by a wand with three reflective markers on it. The maximum acceleration detected when the wand hit the thigh IMU was synchronized with the trajectory of the wand's reflective markers. Data from the OCS were analyzed and the angles were computed in Python. The data from the instrumented knee brace were collected and the angles were computed in Python with a data fusion algorithm based on Valenti et al.

filter applied to the accelerometric and gyroscopic data [START_REF] Valenti | Keeping a good attitude: A quaternion-based orientation filter for IMUs and MARGs[END_REF]. However, data analysis revealed problems with the calibration protocol for the instrumented knee brace due to the dorsal decubitus position of the cadavers. A better alternative estimation was to find the IMU orientation by comparing clusters and IMU movements when there were small flexions in FLEXMAX. This estimation of the IMU orientation was used for all other exercises. Then, some acquisitions were not processed because of the abnormal movement of some of the anatomical marker fixed on the knee brace (MF, LF, MT and LT), which resulted in the selection of only 23 acquisitions in all. Due to motion differences between operators, only 13 acquisitions for pivot shift test were selected among the 23.

Statistical analysis

Statistical analysis was performed on the FLEX100, ROT30 and PS exercises. The variables extracted for reliability were the intraclass correlation coefficient (ICC) and Bland & Altman plots. ICC was considered poor (<0.50), moderate (between 0.5 and 0.75), good (between 0.75 and 0.90), or excellent (>0.90) [START_REF] Cuesta-Vargas | The use of inertial sensors system for human motion analysis[END_REF]. The variables extracted for validity were Pearson's correlation coefficient (σ), the root mean square error (RMSE), and the mean confidence interval (CI). Pearson's coefficient was evaluated as excellent (>0.95), very good (0.85-0.95), good (0.75-0.85), moderate (0.65-0.75) or weak (<0.65) [START_REF] Kayaalp | Validation of a Novel Device for the Knee Monitoring of Orthopaedic Patients[END_REF]. Error (RMSE) <5° was considered excellent, and between 5 and 10° it was considered good [START_REF] Poitras | Validity and reliability ofwearable sensors for joint angle estimation: A systematic review[END_REF]. RMSE was measured without offset on the abduction/adduction (A/A) axis and the internal/external rotation (I/E) axis because of the poor accuracy of the alternative IMU orientation procedure. To deal with STA [START_REF] Teufl | Validity of inertial sensor based 3D joint kinematics of static and dynamic sport and physiotherapy specific movements[END_REF], we extracted two types of variables for the validity: the validity on the standard knee motion and the validity on the knee motion with consideration of STA on the flexion/extension (F/E) axis through a linear regression on FLEXMAX acquisitions.

Results

The ICC (95%) for the kinematic parameters recorded with the instrumented knee brace showed excellent reliability for all exercises. In terms of concurrent validity, Pearson's r coefficient for the instrumented knee brace kinematics compared to the navigation kinematics showed the following for all exercises: very good to excellent correlations for the F/E angles, poor to excellent correlations for I/E rotation, and poor to excellent correlation for A/A (Table 1).

RMSE was less than 5° between systems except for the F/E axis. Table 2, Figure 3, Figure 4 and Figure 5 show the impact of STA. The Bland & Altman plots showed the difference in the range of motion (RoM) between the two systems (Figure 2). 3D kinematics of the knee during the three exercises are depicted in Figure 3, Figure 4, and Figure 5. Rates of flexions for FLEX100 and PS were 78.3 ± 12.4 dps and 55.5 ± 8.9 dps, rates of rotations for ROT30 were 19.4 ± 2.8 dps. Coefficients of variation (CV) were respectively: 15.8%, 16.0%, 14.3%.

Table 1.

ICC and Pearson's correlation coefficient (σ) between instrumented knee brace and navigation.

(CI = confidence interval.)

Discussion

This study compared a new instrumented knee brace based on two IMU sensors aligned to the lateral side of the pelvis and the leg with an OCS to evaluate the accuracy and reliability of the new system during specific tasks. The study showed the excellent reliability of the knee brace and the strong concordance between the two systems. The difference in the RoMs was acceptable. The RMSE values of the OCS and the knee brace for the three angles remained under 5° for all exercises.

The reliability of the IMU system was tested for each exercise and demonstrated ICC>0.918 for all exercises, indicating excellent reliability. Moreover, the results are in line with the literature. Indeed, the ICC of the FLEX100 exercise was similar to the results found by Maderbacher's study, which showed an ICC of 0.92 for knee flexion [START_REF] Maderbacher | Reproducibility of navigation based kinematic analysis of the knee -A cadaveric investigation[END_REF]. The results for I/E rotation with our knee brace were similar to those found by Musahl et al., who reported ICCs from 0.94 to 0.99 for measured rotational knee laxity with a noninvasive system [START_REF] Musahl | Development of a simple device for measurement of rotational knee laxity, Knee Surgery[END_REF]. Other studies of gait analysis using IMU systems have found similar results, with excellent reliability for the flexion angle [START_REF] Cho | Evaluation of validity and reliability of inertial measurement unit-based gait analysis systems[END_REF][START_REF] Kobsar | Reliability of gait analysis using wearable sensors in patients with knee osteoarthritis[END_REF]. These studies are consistent with findings of the excellent reliability of IMU systems for several tasks.

The correlations between the two systems were excellent in the major axis of each exercise.

These results are in line with the Pearson's correlation coefficient of 0.96 reported by Kayaalp et al. and the agreement between 0.82 to 1.0, depending on the task, for knee F/E reported by Lebel et al. [START_REF] Kayaalp | Validation of a Novel Device for the Knee Monitoring of Orthopaedic Patients[END_REF][START_REF] Lebel | Inertial measurement systems for segments and joints kinematics assessment: Towards an understanding of the variations in sensors accuracy[END_REF]. However, the correlations for the minor axes were not as conclusive with

Pearson's correlation coefficients ranging from weak to excellent. Favre et al. brought this trend to light with measurements of IMUs and the knee exoskeleton during a gait test, finding the major axis σ = 1.00 (0.00) in F/E and the minor axes σ = 0.76(0.18) in A/A and σ = 0.85(0.11) in I/E [START_REF] Favre | Functional calibration procedure for 3D knee joint angle description using inertial sensors[END_REF]. Then, additional measurements revealed agreement levels between poor and excellent, depending on the axis and joint, for squats (SQ), single-leg squats (SLS) and countermovement jumps (CMJ). For the knee joint, the minor axes for these exercises showed less agreement, ranging from poor to excellent, compared to the major axis (F/E), which showed excellent agreement [START_REF] Teufl | Validity of inertial sensor based 3D joint kinematics of static and dynamic sport and physiotherapy specific movements[END_REF]. These findings underline the weaker correlation of the two minor axes of the knee and the strong correlation in F/E found in this study, except for the ROT30 exercise due to a major axis in the I/E axis.

The accuracy of the instrumented knee brace as assessed by RMSE was under 5° in the three axes (F/E without STA). Bland & Altman analysis revealed a similar RoM with a mean error under 2° in F/E due to the STA regression used in the computation angles of the instrumented knee brace. These results are similar to the findings of Teufl et al. for SQ, SLS and CMJ with a mean error of 2° [START_REF] Teufl | Validity of inertial sensor based 3D joint kinematics of static and dynamic sport and physiotherapy specific movements[END_REF].

Another study on Sit To Stand (STS) test with visual feedback showed a mean difference ranging from -1.7° to 3.2° in the F/E axis between the IMU system and OCS [START_REF] Bell | Verification of a portable motion tracking system for remote management of physical rehabilitation of the knee[END_REF]. More specifically, Leardini et al. highlighted a mean RoM error of 3.9° for knee flexion between 0 to 95° [START_REF] Leardini | Validation of the angular measurements of a new inertial-measurement-unit based rehabilitation system: Comparison with state-of-the-art gait analysis[END_REF]. However, the RoM error for our ROT30 exercise revealed an underestimation in the I/E rotation axis due to STA. Unfortunately, no significant regression was found to minimize STA on this axis. Indeed, Stagni et al. also observed a simple overestimated flexion in F/E but a more complex artifact for I/E and A/A rotations [START_REF] Stagni | Quantification of soft tissue artefact in motion analysis by combining 3D fluoroscopy and stereophotogrammetry: A study on two subjects[END_REF].

A few limitations of this study should be noted. First, the outcomes were marred by STA, which impacted each axis of rotation, but more specially the I/E and A/A rotations because there was no compensation for STA. More studies are required to identify the best approach to minimize STA and reduce RMSE between the connected knee brace and bone placement [START_REF] Tsai | Quantification of three-dimensional movement of skin markers relative to the underlying bones during functional activities[END_REF][START_REF] Potvin | A practical solution to reduce soft tissue artifact error at the knee using adaptive kinematic constraints[END_REF]. Indeed, the results of the ROT30 exercise contradicted the findings of Moewis et al. (Figure 2; Figure 3), who suggested a linear overestimation of I/E rotation with skin-mounted markers compared to knee fluoroscopy [START_REF] Moewis | Towards understanding knee joint laxity: Errors in non-invasive assessment of joint rotation can be corrected[END_REF][START_REF] Lin | Effects of soft tissue artifacts on differentiating kinematic differences between natural and replaced knee joints during functional activity[END_REF]. Second, the alternative procedure for estimating the IMU orientation with optimization introduced inaccuracy into the measurement of I/E and A/A rotations due to the exercise of flexion used by the optimization algorithm [START_REF] Cooper | Inertial sensor-based knee flexion/extension angle estimation[END_REF][START_REF] Robert-Lachaine | Inertial motion capture validation of 3D knee kinematics at various gait speed on the treadmill with a double-pose calibration[END_REF]. In fact, the I/E and A/A rotations had a small range of motion and were in general less valid due to the projection of the F/E rotation, caused by small errors in IMU orientation estimation during the calibration phases [START_REF] Poitras | Validity and reliability ofwearable sensors for joint angle estimation: A systematic review[END_REF][START_REF] Robert-Lachaine | Inertial motion capture validation of 3D knee kinematics at various gait speed on the treadmill with a double-pose calibration[END_REF][START_REF] Manal | Comparison of surface mounted markers and attachment methods in estimating tibial rotations during walking: an in vivo study[END_REF]. Further studies should elucidate the most appropriate calibration methods to estimate joint kinematics [START_REF] Kayaalp | Validation of a Novel Device for the Knee Monitoring of Orthopaedic Patients[END_REF][START_REF] Cooper | Inertial sensor-based knee flexion/extension angle estimation[END_REF]. Third, even if navigation system allows a motion capture without STA, mechanical properties between cadaveric and living knee may differ due to difference in soft tissue tension, capsule, ligaments and neglecting the influence of muscle contractions [START_REF] Worlicek | The influence of varus and valgus deviation on patellar kinematics in healthy knees: An exploratory cadaver study[END_REF][START_REF] Maderbacher | Influence of Component Rotation in Total Knee Arthroplasty on Tibiofemoral Kinematics-A Cadaveric Investigation[END_REF]. Then, passive kinematics may differ from active kinematics [START_REF] Worlicek | The influence of varus and valgus deviation on patellar kinematics in healthy knees: An exploratory cadaver study[END_REF][START_REF] Maderbacher | Influence of Component Rotation in Total Knee Arthroplasty on Tibiofemoral Kinematics-A Cadaveric Investigation[END_REF]. Not least of all, rates of flexion and rotation for FLEX100, PS and ROT30 highlighted limitations in reliability of movements due to a CV close to or even higher than 15% [START_REF] Schectman | Methods of Clinical Epidemiology[END_REF]. Acquisition performed with a robotic arm would acquire better data with repetitive ranges and rates [START_REF] Du | Effects of Proud Large Osteochondral Plugs on Contact Forces and Knee Kinematics A Robotic Study[END_REF].

However, acquired data performed with clinicians is closer to field data acquisitions.

Conclusion

The present work investigated the reliability and accuracy of a new instrumented knee brace compared to a gold standard navigation system without STA. The new brace showed excellent reliability and sufficient accuracy for clinical interpretations. However, STA induced errors in the axes without artifact consideration algorithm. In conclusion, this study highlights the clinical validity of the instrumented knee brace. Further investigations are necessary to develop new algorithms to minimize STA in the I/E and A/A

Declaration of competing interest

The institutes of the Authors have received funding from the marketing company the

Acknowledgements

The Authors acknowledge Ted Orthopedics for their technical and financial contribution to this study. The study sponsors had no role in the study design and interpretation of data in the writing of the manuscript and in the decision to submit the manuscript for publication.

Figure 1 :

 1 Figure 1: View of the instrumented knee brace, reflective markers and clusters placement on the cadaver's lower limb. The torso was strapped to the table to minimize artifacts on the hip cluster position.

Figure 2 :

 2 Figure 2: Bland & Altman plots of RoM from both systems for FLEX100 (a), ROT30 (b) and PS (c). Each graph represents the mean difference (black line) and the CI 95% of the difference (dashed lines) as recorded by the instrumented knee brace and navigation.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Mean knee joint kinematics over all trials and all subjects with a standard deviation cloud to compare the two systems (Instrumented knee brace in blue, Instrumented knee brace without STA regression in red, Navigation with Vicon in orange) for flexion 100° in F/E (a), I/E (b), A/A(c)

  instrumented knee brace Ted Orthopedics to cover the costs involved in the execution of this study. I, Nicolas Reneaud, declare that Zory Raphaël, Chorin Frédéric, Chavet Pascale, Coyle Thelma, Puech Stephane, Ollivier Matthieu, Chabrand Patrick, Gerus Pauline have no proprietary, financial, professional or other personal interest of any nature of kind of any product, service and/or company that could be construed as influencing the position presented in the manuscript.

Table 2 .

 2

	RMSE of F/E, F/E with STA regression, I/E, and A/A with 95% confidence interval
	Exercise	FLEX100 (°)	ROT30 (°)	PS (°)
	F/E (CI 95%)	8.26 (6.77 -9.77)	4.65 (3.11 -6.18)	6.06 (4.86 -7.27)
	F/E Reg (CI 95%)	2.96 (2.31 -3.61)	1.01 (0.83 -1.20)	2.72 (2.08 -3.35)
	I/E (CI 95%)	2.25 (1.99 -2.52)	3.12 (2.67 -3.58)	3.89 (3.03 -4.75)
	A/A (CI 95%)	3.82 (3.25 -4.38)	1.83 (1.57 -2.10)	3.53 (2.45 -4.61)

  Lebel et al. reported comparable mean RMSE varying between 1.1° and 5.5°, depending on the segment tracked and the task performed[START_REF] Lebel | Inertial measurement systems for segments and joints kinematics assessment: Towards an understanding of the variations in sensors accuracy[END_REF]. The same results were found by Kayaalp et al., with a mean RMSE of 5.17° in F/E during a gait test[START_REF] Kayaalp | Validation of a Novel Device for the Knee Monitoring of Orthopaedic Patients[END_REF]. A systematic review by McGinley et al. confirmed these results and noted that most studies reported error of less than 5° for all gait variables, excluding hip and knee rotation[START_REF] Mcginley | The reliability of three-dimensional kinematic gait measurements: A systematic review[END_REF]. The study by Robert-

	Lachaine et al. was in line with these results, with RMSE between 3° and 6° for knee rotations
	in the F/E, I/E, and A/A axes [51]. These authors thus agreed on the clinical relevance of IMU
	systems with RMSE values under 5°.