Nourddine Azzaoui

Member, IEEE Laurent Clavier

UWB channel modeling for objects evolving in impulsive environnements

We consider channel modeling issues in the context where communicating objects are evolving in impulsive environments. It was shown recently that α-stable random processes are attractive solution for representing the ultra wide band communication channel in relatively large spatial areas. In this paper, we consider the α-stable channel modeling in an evolutionary context where the model features depend on spatial locations. We introduce a methodological approach consisting of two parametric and non parametric components: the latter can be seen as black box model to describe the spatial evolution and it can be learned from historical observations of the transfer function. The other component concerns the frequency dependence and has an auto-regressive structure.

INTRODUCTION

One important difficulty of statistical channel modeling (especially when ultra wide band (UWB) is considered) resides in its ability to represent the time and environment evolutions. This channel variability is obvious when it comes to model mobiles and rapidly changing environments. In order to make realistic simulations, it is necessary to adapt the existing models to such situations. Indeed, statistical models, even those taking into account rare events (α-stable models), are not sufficient to describe the complexity of the channel behavior in all circumstances. Another challenge for developing real world applications is the fact that the specification of a such model needs a large measurement campaign and usually takes a lot of time to estimate the model parameters. For the next generation communication systems, (as WSN, inter vehicles, ...), time evolution can be slow or fast and nodes will change their position and interact with different communication infrastructures. They must be able to learn their environments autonomously, especially channel models of the medium in which they evolve. Due to the models complexity and possible fast evolution, this task is quasi impossible for the end user, especially low complexity sensor nodes.

The idea of this work is to give a general model that can adapt to the rapid change of the environment and learned from a limited number of measurements. For this purpose our proposed model must, from one hand, describe the space dependence of the channel transfer function and, on the other hand, the model variation with frequency. From this last point of view many works described this frequency dependence as Nourddine Azzaoui is with the Mathematics Laboratory UMR-CNRS 6620 from of Blaise Pascal university, Aubière (France).

Laurent Clavier is with IEMN (Institute of Electronics, Microelectronics and Nanotechnology, UMR CNRS 8520), IRCICA (Research Institute on software and hardware components of the future for information and communication, FR CNRS 3024) and TELECOM Institute, TELECOM Lille 1, Lille, France. This work was supported by the ERDF (European Regional Development Fund) and by the Nord-Pas-De-Calais Region (France). parametric model: we cite for instance Ghassemzadeh [START_REF] Ghassemzadeh | Measurement and modeling of an ultra-wide bandwidth indoor channel[END_REF] who uses a second-order autoregressive model AR(2) for frequency response generation of the UWB indoor channel. The work in this paper is an on going contribution, it presents the theory behind the model and solutions for space and frequency channel evolution in section I and some mathematical tools to estimate the parameters (section II) and validate the model (section III). Further investigations is needed to validate this approach. We also assume that the channel is accurately represented by an α-stable model as we proposed in [START_REF] Azzaoui | Statistical channel model based on αstable random processes and application to the 60 GHz ultra wide band channel[END_REF].

I. THE EVOLUTIONARY CHANNEL MODEL

The main idea of the evolutionary model comes from the fact that transfer functions measured at two contiguous places will be strongly dependent; it will be weakly dependent when the distance increases between two measurements. For the mathematical formulation let us denote H x the transfer function at a position x. We postulate that the dependence will have a memory of length p. This space dependence can be formulated using conditional expectations as follows:

H x+1 = f (H x , . . . , H x-p ) + Υ, (1)
where f is an unknown function not depending on the location and Υ = H x+1 -E (H x+1 | H x , . . . , H x-p ). Equation (1) can be seen as a black-box model which does not support any a priori about the environment. One of its main advantages is the fact that it can be learned (estimated non parametrically using kernel techniques) from historical observations collected by a given node. In order to understand the frequency dependence we present the model for observed transfer functions at w 1 , . . . , w n , as the following:

pos p      H p+1 (w 1 ) = f (H p (w 1 ), . . . , H 1 (w 1 )) + Υ p (w 1 ) . . . . . . H p+1 (w n ) = f (H p (w n ), . . . , H 1 (w n )) + Υ p (w n ) . . . . . . . . . pos x      H x+1 (w 1 ) = f (H x (w 1 ), . . . , H x-p (w 1 )) + Υ x (w 1 ) . . . . . . H x+1 (w n ) = f (H x (w n ), , . . . , H x-p (w n )) + Υ x (w n )
Index p denotes the length of the dependence, consequently, the number of samples that nodes has to store in order to start its estimation. For different positions x and x , we suppose that vectors Υ x and Υ x are pairwise independent. Since the function f has captured all the intrinsic spatial influence, it makes more sense to suppose that Υ depends only on frequencies. Inspired from works in [START_REF] Ghassemzadeh | Measurement and modeling of an ultra-wide bandwidth indoor channel[END_REF], we suppose that for a fixed position, the internal frequency dependence of Υ may be described by an autoregressive scheme.

Υx(wi) = q k=1 θ k Υx(w i-k ) + εx,i (2)
In order to take into account the impulsive nature and the α-stable model [START_REF] Azzaoui | Statistical channel model based on αstable random processes and application to the 60 GHz ultra wide band channel[END_REF], we suppose that errors ε x,i 's are i.i.d. symmetric α-stable white noises.

The main difficulty in the black-box model ( 1) arises when we try to estimate the non parametric function f which will be seriously affected when p becomes large; this is known as the curse of dimensionality. In order to overcome such an inconvenience we will suppose that, instead of (H x , . . . , H x-p ), the transfer function H x+1 depends on a linear combination of this vector. This reduce f to a real function and hence its estimation is facilitated. This will lead to the following:

Hx+1(wi) = f ( p k=0 η k H x-k (ωi)) + q k=1 θ k Υx(w i-k ) + εx,i (3)
To resume, f represents the space dependence, the coefficient θ k the frequency linear dependence and ε x,i the unpredictable variations. This model is very general and will allow the use of powerful mathematical tools known in literature as semi parametric single index partially linear model. They have been largely studied in literature we cite among others [START_REF] Carroll | Generalized partially linear single-index models[END_REF], [START_REF] Gao | Statistical inference in single-index and partially nonlinear models[END_REF], [START_REF] Xia | On extended partially linear single-index models[END_REF], [START_REF] Xia | A goodness-of-fit test for singleindex models[END_REF]... For compatibility with the semi-parametric formalism we rearrange the observed transfer functions as follows: we begin by taking N = nm where n is the number of positions and m is the number of observed frequencies. For every x = 1, . . . , n -1, for every i = 1, . . . , m we take s = (x -1)m + i and we denote:

Y s = H x+1 (w i ), U s = [H x (w i ), . . . , H x-p (w i )] † , V s = [Υ x (w i-1 ), . . . , Υ x (w i-q ) ] † , X s = [U † s , V † s ] † , ε s = ε x,i
where † states for matrix transpose. This will lead to the conventional semi parametric single index model notations:

Ys = X † s θ + f (X † s η) + εs (4)
where θ = (0 p+1 , θ 1 , . . . , θ q ) and η = (η 0 , . . . , η p , 0 q ) and 0 p is zeroes vector in R p . Motivation for using (4) for independent data analysis can be found in [START_REF] Xia | On extended partially linear single-index models[END_REF], [START_REF] Xia | A goodness-of-fit test for singleindex models[END_REF], [START_REF] Gao | Nonlinear Time Series Semiparametric and Nonparametric Methods[END_REF]. Many indications suggest the use of semi-parametric models (4); for example two transfer functions are presented in Fig. 1. Figures Fig. 2 and Fig. 3 illustrate the space and frequency dependence. The autoregressive scheme was also noticed in [START_REF] Ghassemzadeh | Measurement and modeling of an ultra-wide bandwidth indoor channel[END_REF] and confort our proposal.

II. STATISTICAL ESTIMATION AND MODEL SPECIFICATION FOR THE EVOLUTIONARY MODEL

In this paper, we will address the estimation problem in the α-stable case and introduce consistent estimators of unknown model (4) features; the mode of convergence and the quality of the estimators will be detailed in further works.

A. Estimation in semi parametric single index models

Consider a semi-parametric single-index model of the form:

Y s = X † s θ + f (X † s η) + ε s , s = 1, 2, . . . N, (5)
where both θ and η are vectors of unknown parameters. The real function f (•) is unknown defined on R and is supposed to be twice differentiable. The random variables {ε s } are a sequence of errors with E[ε s | X s ] = 0. A huge amount of works about such models have been introduced in literature especially when the errors are random variables with finite second order moments; for more details see for instance [START_REF] Xia | On extended partially linear single-index models[END_REF] and [START_REF] Gao | Nonlinear Time Series Semiparametric and Nonparametric Methods[END_REF]. However no works have been done in the case of infinite variance errors. In this paper we suppose that the ε s 's are i.i.d. α-stable centered variables with a fixed scale parameter σ.

In order to estimate the unknown parameters η, θ and the function f involved in (5), we introduce the following notation:

f 1η (u) = E[Ys | X † s η = u], (6)
f 2η (u) = E[Xs | X † s η = u], (7)
This decomposition is inspired from the conditional decomposition

f (u) = E[Y s -X † s θ | X † s η = u]
which leads easily to the formula:

f (u) = f 1η (u) -(f 2η (u)) † θ (8)
Let K be a kernel function i.e. a non negative even probability density (for example the gaussian kernel K 2 ). From the conditional expectation given in [START_REF] Xia | A goodness-of-fit test for singleindex models[END_REF] we propose a kernel estimation inspired from the Nadarya-Watson estimation techniques. We first estimate f 1η (•) by:

(x) = 1 √ 2π e -x 2
f1η (u) = N s=1 K h (X † s η -u)Ys N s=1 K h (X † s η -u)
where K h (•) = K( • h ) and h is a bandwidth parameter. Similarly we estimate the function f 2η (•) by:

f2η (u) = N s=1 K h (X † s η -u)Xs N s=1 K h (X † s η -u)
For a complete estimation of the function f we will need an estimate of the parameters η and θ. For this purpose we use the following notations:

Y sη = Ys -f1η (X † s η), Xsη = Xs -f2η (X † s η),
Let us consider the least-squares sum :

SN (θ, η; h) = N s=1 (Y sη -X † sη θ) 2
The estimation idea consist in minimizing S N (θ, η, h) over (θ, η, h). We first remark that, for a fixed (η, h), the least squares estimator of θ can be deduced using classical linear regression techniques, it is given by:

θ(η, h) = N s=1 XsηX † sη + N s=1 XsηY sη , (9)
where ( . ) + denotes matrix pseudo-inverse. We then estimate (η, h) by (η, ĥ) through minimizing,

ŜN (η, h) = N s=1 (Y sη -X † sη θ(η, h)) 2 (10)
Inspired from the equation ( 8), we propose the nonparametric estimator of f (•) by using estimates of f 1η and f 2η :

f (u) = f1η (u) -( f2η (u)) † θ(η, ĥ). (11)
On the other hand, when the errors scale parameter σ is unknown, it can be estimated using the fractional lower moments technique as follows:

σ = Cα(ρ) 1 N N s=1 Y sη -X † sη θ(η, ĥ) ρ 1 ρ (12)
for every 1 < ρ < α and C α (ρ) is an universal constant depending only on α and ρ, it is given by:

C α (ρ) = α √ π Γ( -ρ 2 ) 2 ρ+1 Γ( 1+ρ 2 ) Γ( -ρ α ) 1 ρ
The mode of convergence and the consistency of these estimators was studied in the second order case, the reader find a detailed overview in literature; see for instance [START_REF] Gao | Nonlinear Time Series Semiparametric and Nonparametric Methods[END_REF]. In this last work a central limit theorem type result was established for all estimators. We believe that similar results may be proven for the α-stable process and will be similar to the generalized central limit theorem context.

III. MODELS ADEQUACY AND HYPOTHESIS TESTING

In this section we present statistical techniques to test the adequacy of semi parametric models of type [START_REF] Gao | Statistical inference in single-index and partially nonlinear models[END_REF]. Recently, semiparametric approach has been used for model specification tests in the case of finite variance processes. We will proceed by analogy to second order works, we focus on partially linear or single-index scheme against a general non-parametric form. We concentrate on parametric specification testing of the conditional mean function defined for u ∈ R p+q+1 by:

m(u) = E[Y s | X s = u)]

A. Testing for single-index regression

In this subsection we focus on the particular case of (4) without the linear regression component. The test purpose is to see if transfer functions evolutions can be reduced to the space evolution, not considering the frequency dependency (H 0 ). Somehow it can be seen as a test for spatial stationarity of transfer functions. We thus look at testing the null hypothesis of single index modeling against a class of non parametric functions:

(H 0 ) : m(x) = f (x † η), (H 1 ) : m(x) = f (x † η) + ∆(x) for all x ∈ R p+1 ,
where f (•) is an unknown function on R, η is a vector of unknown parameters and ∆ is any regular function defined on R p+1 . Under the null hypothesis (H 0 ) we have given techniques to estimate η in section II. In the single index particular case the estimate of f (•) will be simplified and is given by:

f (X † s η) = N t=1 YtK( (Xs -Xt) † η h ) N u=1 K( (Xs -Xu) † η h ) ,
with K(•) being a kernel function defined on R. Consequently, from [START_REF] Stute | Nonparametric checks for single-index models[END_REF] the parameter η is then estimated by minimizing:

η = arg min (η,h) N s=1 (Ys -f (X † s η)) 2
As we have supposed that the errors ε s 's are α-stable i.i.d. centered random variables, then by the generalized central limit theorem it is more convenient to approximate the estimated errors Ŷs = Y s -f (X † s η) by a centered symmetric α-stable distribution with estimated scale parameter σ:

σ = C α (ρ) 1 N N s=1 Ŷs ρ 1 ρ
By analogy with the second order case and for compatibility with infinite variance we suggest the test statistic:

L = N s=1 | Ŷs | σ
Many improvements are still needed in the case of stable variables, including the asymptotic distribution and robustness of the statistic L. For the hypotheses testing in the case of second order semi parametric models one can find a rich literature in [START_REF] Fan | Consistent model specification tests: omitted variables and semiparametric functional forms[END_REF], [START_REF] Gao | Statistical inference in single-index and partially nonlinear models[END_REF], [START_REF] Stute | Nonparametric checks for single-index models[END_REF] and references within.

B. Testing for partially linear single-index model

It is a natural extension of the single index model to situations where the linear regression component may be suitable for modeling the inherent studied phenomenon. The following test evaluates if the dependency on space and frequency is sufficient (H 0 ) or if it would be more accurate to take some more phenomena into account, for instance some non linear dependency in frequency.

(H 0 ) : m(x) = x † θ + f (x † η) (H 1 ) : m(x) = x † θ + f (x † η) + ∆(x),
for all x ∈ R p+q+1 . This test problem has been studied for the second order case by [START_REF] Xia | On extended partially linear single-index models[END_REF], [START_REF] Xia | A goodness-of-fit test for singleindex models[END_REF]. With the same reasoning as in the single index model and by using the estimation techniques under the hypothesis (H 0 ) presented in section II, we propose the test statistic:

L = N s=1 | Y s | σ
where Y s = Y s -X † s θf (X † s η), the quantities θ, η, f (•) and σ are respectively the consistent estimators given in ( 9), ( 10), (11) and (12). One of the main issues to evaluate the models compatibility is to examine the power and the size tests : α = P(L > l α | such that H 0 "holds" )

β = P(L > l α | such that H 1 "holds" )
where l α will be tabulated from the distribution of L. It should be noted that the estimator presented here must be extensively studied for asymptotic properties and robustness.

CONCLUSION AND PERSPECTIVES

In this paper we have introduced a methodological approach that adapts semi-parametric techniques to the context of spatial ultra wide band channel modeling. Our approach may be seen as a general-purpose framework for the modeling of the spatial evolution of the channel transfer functions described by an αstable random process. Extensive simulations on real data must be performed to evaluate the adequacy of the proposed model in the context of indoor and outdoor configurations and this paper only presents the general framework.

Fig. 1 .Fig. 2 .

 12 Fig.1. Examples of measured transfer functions between 57-59 GHz and down converted between 1-3 GHz at two different locations. More details on the measurement setup can be found in[START_REF] Azzaoui | Statistical channel model based on αstable random processes and application to the 60 GHz ultra wide band channel[END_REF],[START_REF] Fryziel | Path loss model of the 60 GHz radio channel[END_REF] .

Fig. 3 .

 3 Fig. 3. Frequency dependence structure in observed transfer functions. The frequency step is 5 MHz.


Abdelkader Miraoui
email: abdelkader.miraoui@utt.fr

Hichem Snoussi

Jacques Duchêne

Nourddine Azzaoui

On the Detection of Elderly Equilibrium Degradation Using Multivariate-EMD

The aim of this paper is to provide a new methodology for the detection of an increased risk of falling in community-dwelling elderly. A new extended method of the empirical mode decomposition (EMD) called multivariate-EMD is employed in the proposed solution. This method will be mainly used to analyze the stabilogram center of pressure (COP) time series. In this paper, we describe also the remote non-invasive assessment method, which is suitable for static and dynamic balance. Balance was assessed using a miniature force plate, while gait was assessed using wireless sensors placed in a corridor of the home. The experimental results show the effectiveness of this indicator to identify the differences in standing posture between different groups of population.

I. INTRODUCTION

Due to the resulting loss of the elderly autonomy, the deterioration of their balance is a major problem in biomedical. Thus, several risk factors have been identified to minimize the effects of the falls. The muscular weakness is recognized as the most commonly cited factor. Recent analysis [START_REF] Prieto | Measures of postural steadiness: differences between healthy young and[END_REF]- [START_REF] Brauer | A prospective study of laboratory and clinical measures of postural stability to predict community-dwelling fallers[END_REF] show that the economic cost of the falls is proportional to the elderly population. In France, for example, the annual number of deaths caused by the falls is estimated to be more than 9.000, with a resultant cost over than two billion euros [START_REF]Prévention des chutes chez les personnes âgées : des conseils pour préserver l'autonomie et le bienêtre de nos aînés[END_REF].

To capture clinically relevant episodes accurately and reliably, in a pervasive healthcare monitoring system, multiple sensors are employed to measure, both physiological and contextual informations. The postural equilibrium is maintained by reacting information from all the sensory systems, including vestibular, visual, and proprioception systems.

To evaluate postural sway, a force plate is used which integrates wireless sensors recording data captured from the body mouvements. The force plate measures the displacement of the COP, which represents the location of the resultant force exerted on the surface of a force plateform. The COP, which can be used as a measure of postural stability, is defined in the horizontal plane in both anteroposterior (AP) and mediolateral (ML) directions [START_REF] Prieto | Measures of postural steadiness: differences between healthy young and[END_REF]. The representation of the COP time series, in the AP and ML directions, is known as the stabilogram (see Fig. 1).

Authors in [START_REF] Carroll | Nonstationary properties of postural sway[END_REF] and [START_REF] Harris | Standing postural stability assessment and signal stationarity in children with cerebral palsy[END_REF] have proved that the stabilogram signals are nonstationary. Hence, the use of these curves presents a lack of ability to decide whether the person has a risk of falling or not. Therefore, we need a more precise and powerful method to explore these signals and then extract the relevant parameter. In this paper, in order to extract informations from the stabilogram and to characterize the postural stability, we propose the use of a new version of empirical mode decomposition of bivariate signals. This is the stabilogram signal, which may be seen as a multi-component signal [7]. Fig. 1 depicts an example of real data scenario.

The proposed method ensures the decomposition of the signals into elementary rotating functions, oscillating functions and tendencies. The technique is based on recursive extraction of non planar and rapidly rotating functions, and then a recursive extraction of rapidly oscillating functions as in case of univariate-EMD. Our concept is based on a theoretical decomposition and is confirmed on simulated signals (see Fig. 2).

We propose to exploit the stabilogram to extract a relevant indicator that allows diagnosing a risk of falling and monitoring of the evolution of imbalance quality in the elderly.

The remaining of the paper is organized as follows. Firstly, the experimental protocol, mv-EMD, and the data processing are presented in Section II. Section III summarizes the main result, obtained by applying our proposed method on the stabilogram. Finally, Section IV concludes the paper. In this section, we give a detailed description of the proposed method and the materials employed for signal acquisition and data analysis.

A. Population Groups and data Acquisition

During all manipulations, we assume that all measurements are taken from two classes of people: the first class is composed of ten healthy control subjects (three males and seven females) while the second includes ten healthy elderly subjects (four males and six females).

Control subjects mean age, height, and weight are 33.3±7.4 years, 168.0±6.5 cm, and 65.7±17.6 kg, respectively. Elderly subjects mean age, height and weight are 80.5 ± 4.7 years, 165.6 ± 7.0 cm and 71.9 ± 9.9 kg, respectively.

We notice here, that the data related to center of pressure (COP) are obtained from a Bertec 4060-08 force plate (Bertec Corporation, Columbus, Ohio, USA). The initial COP signals are recorded for several non-normalized values of the center of force plate and by subtracting the mean value. The data is sampled at the rate of 100 Hz and then filtered using an 8thorder lowpass Butterworth filter with a cutoff frequency of 10 Hz. Then the data are recorded and treated using ProTags and Labview (National Instruments Corporation, Austin, Tex, USA). Finally, all subsequent calculations are performed using MATLAB.

B. Experimental protocol

Subjects are tested barefoot or wearing socks. The test starts with subjects standing upright with their arms by their sides in front of the force-plate while looking at a 10 cm cross fixed on the wall in 2 meters in front of them. Upon verbal instructions, subjects stepped onto the force plate and they are not required to use a preordained foot position. Data recording lasted 15 seconds, during which subjects maintain an upright posture.

C. Data processing 1) Overview of the mv-EMD:

The EMD concept was first introduced by Huang et al. [7]. It is based on the simple assumption that any data may have many different coexisting modes of oscillations one superimposing on the others. It is a signal processing decomposition technique that decomposes each signal into waveforms modulated in both amplitude and frequency by extracting all modes that the signal contains [7].

The decomposition is an intuitive and adaptive signaldependent and does not require any conditions about the stationarity and linearity of the signal. The waveforms extracted by EMD are named IMFs (Intrinsic Mode Functions). Each IMF is symmetric, assumed to yield a meaningful local frequency, while different IMFs do not exhibit the same frequency at the same time. In other words, each IMF satisfies the two following constraints:

• the number of extrema and the number of zero crossings are identical or differ at most by one;

• the mean value between the upper and the lower envelope is equal to zero at any time. The difference between the original signal and the IMF time series is the residual. The first IMF component is obtained by a sifting process. This procedure is then applied on the residual in order to extract the second IMF, and so forth. Thus all the IMFs are iteratively extracted. The nonstationary signal ( ) is then represented as a linear sum of IMFs and the residual component. After the convergence of this algorithm the signal ( ) can then be written as:

( ) = ∑ ( ) + ( ), (1)
where ( ) is a residual tendency and ( ) are the intrinsic modes IMFs. (IMF : intrinsic mode function) It relies heavily on the calculation of the extrema that do not have any meaning in such multidimensional signals.

The originality of our method lies in the fact that it drops the step of sifting and consequently its drawbacks. In contrast to the classical EMD, our proposed method uses the convexity changes of the signal instead of its extrema. These new adaptations are easy to deal with, since they have no natural significance for vectorial or complex functions. Some research works related to the applications of EMD, in the case of vectorial and complex signals, are provided by authors in [8],

[9] respectively .

Our proposition is to be able to decompose any signal into elementary bivariate IMFs which is a bivariate function characterized having the following forms :

• Rotating IMFs: are non-planar curves which are turning around the time axis (they have a spiral shape). In addition, they are "locally symmetric" as explained bellow.

For example (sin( ), cos( )) is a non planar rotating IMF.

• Oscillating IMFs: are planar curves having many changes of the convexity and have a "local symmetry" as well. In their containing plan, they can be seen as a univariate IMFs. For example, the function (sin( ), 2 sin( )) is a simple planar oscillating IMF.

• Tendencies: are planar or non planar curves which have no inflexion. What will be considered as tendency is any function which has no change in the sense of its convexity. Somehow, a tendency may be seen as a piece of an IMF that one cannot observe entirely because of the relative short time of observations. In order to decompose the signal into elementary IMFs, the procedure consists on identifying the inflexion points by searching the points for which the convexity of the curve changes. These inflexions give an idea about the apparent rotation and oscillatory nature of the studied signal.

2) Model motivation and physical explanation: The balance of a person, moving or in a static state, can be described by the fundamental rule of dynamics, illustrated by a differential equation of Langevin type:

⃗ ( ) ( ) = -( ) . ⃗ ( ) + ⃗ ( ), (2)
where ⃗ = ( ˙ , ˙ ) is the velocity vector, and ⃗ ( ) is the resultant of forces applied by the muscles to maintain the balance. The term ( )⃗ ( ) represents the resultant forces responsible for the imbalance. In [7], it was shown that the empirical mode decomposition is more suited to this kind of system (2) in its one-dimensional version. It breaks down the movement into elementary displacement which can be described by intrinsic components (IMFs). In this paper we use the bi-dimensional version of the EMD introduced in [10] which will consider the mutual link between and and their combined influence on the balance. We are aware that the EMD will not necessary find exactly the real elementary behaviors that gave rise to the movement. However, it will give exploitable intrinsic mode functions summarized through their instantaneous amplitudes and frequencies. These instantaneous components will allow the construction of an accurate indicator permitting to point out the state of balance degradation of a given patient.

3) An indicator for precarious imbalance detection:

The idea is based on a physiological explanation of the equilibrium. Firstly, in order to maintain equilibrium, muscles must deploy more energy proportionally to the gap between the COP and the center of gravity (COG). This means that equilibrium becomes more precarious as well as this gap is larger. This distance between the COP and COG will manifest as an amplitude in the stabilogram. On the other hand, the precarious state of the equilibrium worsens when the patient spends more time in this situation of imbalance. This fact will be seen as a the time period to return to the state where the COP coincides with the COG. The idea of the indicator we introduce here consists in exploiting these periods and amplitudes. These last quantities cannot be directly extracted from the stabilogram but only from its corresponding IMFs. For this reason, we begin first by decomposing the stabilogram using the mv-EMD and then extract instantaneous amplitudes and frequencies. The indicator that allows to quantify the risk of falling is introduced as a deviation from the equilibrium amplified by the time spent in this imbalance state.

Indeed, for each IMF extracted from the stabilogram and each simple oscillation, the parameter we propose is given by = . , where is the amplitude of a given local oscillation and its corresponding period. The amplitude , weighted by the period , quantifies the influence of the deviation from equilibrium which is amplified by the duration of the imbalance. In practice, especially when it comes to rapid oscillation IMFs, it is unnecessary to process the data for all periods but rather at fixed laps of time that we choose. The indicator has been extracted using the forthcoming rules:

• Data acquisition and extraction of the stabilogram.

• Signal decomposition using the mv-EMD and extraction of all the IMFs and trends (Suppose we obtain IMFs).

• The indicator , = 1 . . . , is calculated for all the IMFs and for every local oscillation.

• The interval of 10 seconds is fractionated into small intervals. Then, we calculate max the maximum value of the indicator over an interval of 2 seconds corresponding to the th IMF.

• The final indicator, represented in figure (Fig. 9), is equal to = ∑ =1 max .

• We repeat the same process for each subject and for each category.

III. SIMULATION RESULTS

The idea of control (or visual performance) is inspired from technique of management and quality control which are classically used in manufacturing processes and systems reliability. This is achieved by establishing a control chart (Fig. 9) which is composed of a curve representing the evolution of the indicator for ten seconds over three regions: a normal, alert and risk zones. The bounds of these regions are obtained from a large control population. This control chart reflects the ability of an older person to maintain his balance and eventually detect the risk of fall.

A. Signal decomposition

For the implementation of our technique, we illustrate its ability to separate our data. For this purpose, we propose to study an interval of 2 seconds of the real signal. In Fig. 3 to 8, we give the curve shape of an example and its components respectively. We see clearly that the technique introduced in this paper is able to decompose the signal into IMFs and separate rapidly rotating and oscillating components.

B. Extraction of the indicator

Taking curve corresponding to a representative population of young controls as a reference and a large number of autonomous and non autonomous older persons; we have established the chart given in Fig. 9. In this chart we

AP deplacements (cm) ML deplacements (cm)

Fig. 4. the shape of the extracted first IMF distinguish three zones delimited by two bounds: the lower limit corresponds to the the maximum value of the indicator corresponding to the worst situation of young subjects. The dashed line corresponding to the limit of the risk zone is obtained from the worst case corresponding to autonomous but older persons. This way we have sketched a reference graph that shows the variation bounds of the indicator versus time. Fig. 9, shows the temporal evolution of the indicator. It clearly shows that the curve corresponding to the young controls is less chaotic than the elderly. This is quite logical because of the limited autonomy of the latter. Furthermore, the difference of behavior between elderly and autonomous persons is clearly highlighted by our indicator. The detection of a fall risk is reflected simply by the fact of finding a points outside these limits. In this paper we give a simple indicator which can put in evidence the state of equilibrium of vulnerable or older persons. This is based on a novel technique of signal processing using multivariate EMD. This signal decomposition method is applied to the stabilogram whose goal is to give a classical tool that allows diagnosing a risk of falling. Besides its simplicity, the proposed indicator alow monitoring the evolution of balance quality in the elderly. This is important since it can be used to detect or prevent some pathology linked to equilibrium in elderly populations. Another advantage is that this indicator can be adapted to a real-time tracking of the state of the equilibrium and can be incorporated into prosthesis systems for real time monitoring.

Fig. 1 .

 1 Fig. 1. Time evolution of ML and AP in young people and elderly

Fig. 2 .

 2 Fig. 2. Schematic diagram of a remote monitoring system for balance assessment

Fig. 3 .

 3 Fig. 3. stabilogram signal lasting 2 seconds before the decomposition

Fig. 5 .Fig. 6 .

 56 Fig. 5. the shape of the extracted second IMF

Fig. 7 .Fig. 8 .

 78 Fig. 7. the shape of the extracted fourth IMF


Simon Boyé
email: boye@labri.fr

Pascal Barla

Gaël Guennebaud

A Vectorial Solver for Free-form Vector Gradients

Keywords: Vector graphics, diffusion curves, finite elements

: Free-form vector gradients: Our vectorial solver permits to create complex vector gradients (a). Our method is based on a triangular representation (b) that is output-insensitive and thus works with arbitrary high image resolutions. Our solver does not need to be updated for a variety of operations: (c) instancing, layering and deformation; (d) texture mapping; and even (e) environment mapping.

Introduction

Vector graphics techniques have been tremendously improved in recent years, both through research and software products. Two different types of applications have motivated this burst of innovation: image vectorization and free-form drawing. Vectorization methods usually require a dense set of control points to accurately represent all the details of an input image. With drawing tools, this is quite the contrary: a sparse set of vector primitives with rich expressivity is preferred, since it provides for more direct editing and creation. In this paper, we are primarily interested in free-form vector drawing.

One of the first vector primitives that made tractable the creation of complex color gradients is the Gradient Mesh [START_REF] Sun | Image vectorization using optimized gradient meshes[END_REF]] (GM). It takes inspiration from geometric modeling, and lets artists specify color gradients by assigning color values and derivatives to the vertices of a 2D control mesh. With increasing mesh resolution, this technique permits to represent very complex color gradients. On the down side, the denser the mesh, the more difficult it becomes to edit it (Section 2.1).

Recent work has shown that complex color gradients can advantageously be produced through the harmonic or bi-harmonic interpolation of a few spline curves and sample points defining colors [START_REF] Orzan | Diffusion curves: A vector representation for smooth-shaded images[END_REF][START_REF] Finch | Freeform vector graphics with controlled thin-plate splines[END_REF]. Such Diffusion Curves (DC) based representations have also been successfully used for the design of 3D images, height fields or normal maps (Section 2.2).

Although DC-based representations provide an appealing alternative to GMs, computing such images involves a global differential equation which is rather challenging when seeking for real-time performance, high-quality, and robustness. Current methods strive to rasterize DCs directly in the regular pixel grid of the destination image using finite differences. As explained in more details in Section 2.3, such an approach is subject to aliasing when multiple curves pass through a single pixel, and zooming-in can only be approximated via expensive multiple evaluation passes thus breaking the infinite resolution property of the vector representation. Moreover, the expensive differential equation has to be fully solved not only when the curves are edited, but also for simple operations such as translation, rotation and scaling of the entire drawing. This is a major practical limitation, since it makes DCs not tractable for applications that require instancing, deformation or 3D mapping.

In this paper we address all these limitations by introducing a novel vectorial solver for DC images, i.e., a solver that outputs a vectorial representation. The key idea is to dynamically convert a DC image into a high-order mesh-based representation that is automatically adapted to the complexity of the input. To this end, we propose a robust solution that builds on constrained Delaunay triangulations and Finite Element Methods (FEM). Contrary to mesh-based approaches, this intermediate representation is hidden from the user, and updated only when DCs are edited; hence the approach takes the best of both worlds. A direct benefit is that a variety of common operations do not require any update of the solution, as illustrated in Figure 1a-c. Moreover, since the solution is vectorial, it may be accessed randomly, permitting novel uses of free-form vector graphics such as texture or environment mapping (Figure 1d-e).

In addition, we introduce the notion of asymmetric transmission that permits to progressively and precisely limit the influence of some primitives, thus allowing for the first time the mix of global and local diffusion. We present this extension in Section 3, together with a unified taxonomy of curve constraints that identifies novel types of artistic controls. Our vectorial solver and intermediate representation are detailed in Section 4, and we show in Section 5 a variety of applications leveraging their benefits. Finally, we discuss limitations as well as further extensions in Section 6.

Previous works

The focus of this paper is on free-form vector drawing, which among other requirements demands sparse controls, interactive feedback and fast rendering. We review methods that target these goals, and may be used either for creating vector graphics from scratch, or for editing images converted to vector representations.

Mesh-based techniques

One of the earliest vectorization methods is the ArDeco system [Lecot and Lévy 2006]: it uses a mesh-based representation to convert bitmap images into linear, radial or quadratic gradients. In the same vein, GMs consist of cubic color patches arranged into a quad mesh structure [START_REF] Sun | Image vectorization using optimized gradient meshes[END_REF]]. More recently, [START_REF] Liao | A subdivision-based representation for vector image editing[END_REF] have proposed to use triangular subdivisions surfaces to relax the constraints imposed by the regular structure of GMs. However, all these representations are rather difficult to create from scratch or even to manipulate since they expose to the user a lot of control points, and require to explicitly deal with the underlying mesh structure. This is why such methods have been often confined to automatic conversion.

Our approach also makes use of a triangular mesh, but it serves as an intermediate representation that is automatically generated and hidden from users.

Diffusion Curves: primitives

DCs have been introduced to provide a sparse alternative to GMs. Among other benefits, their number may be adapted to the amount of details one desires to produce, they do not have to be closed, and position and color control points do not have to be co-located. They have been defined as Bézier curves with color controls on each side of the curve [START_REF] Orzan | Diffusion curves: A vector representation for smooth-shaded images[END_REF]. Colors are interpolated along each curve, and then propagated to the rest of the image by solving a Laplace equation. Blurring may be applied in post-process to smooth transitions across curves. Curves and color/blur control points may be set by hand, or captured from an image using for instance the optimized solution of [START_REF] Jeschke | Estimating color and texture parameters for vector graphics[END_REF].

The artistic control of DCs has been improved by several extensions such as blockers and directional diffusion [START_REF] Bezerra | Diffusion constraints for vector graphics[END_REF][START_REF] Bowers | A ray tracing approach to diffusion curves[END_REF]. [START_REF] Finch | Freeform vector graphics with controlled thin-plate splines[END_REF] proposed to solve the higher order Bi-Laplace equation to avoid "tent-like" gradients and produce a smoother color interpolation.

DCs have been used for the creation of other types of gradients: Winnemoeller et al. [2009] propagate normal coordinates from silhouettes and creases, normalizing them in post-process; [START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF] create height fields by diffusing height from curves that represent ridges and cliffs; [START_REF] Jeschke | Estimating color and texture parameters for vector graphics[END_REF] use diffusion to control the parameters of procedural textures. DCs have also been extended to 3D by Jeschke et al. [2009b] to render high quality surface details on 3D objects and by [START_REF] Takayama | Volumetric modeling with diffusion surfaces[END_REF] to create volumetric models with 3D color gradients.

In our approach, we solve the Bi-Laplace equation, expanding artistic controls with novel constraints; in particular, we provide local control over vector gradients thanks to novel transmission curves, and handle complex deformations like environment mapping.

Diffusion Curves: solvers

Representations based on DCs require specific methods to solve their associated partial differential equation. Interactive perfor- mance is obtained by adopting multi-scale approaches inspired by the multigrid algorithm [START_REF] Briggs | A multigrid tutorial[END_REF]] that applies Jacobi iterations in a coarse-to-fine scheme. Despite GPU implementations of the multigrid, the convergence may still be slow to generate images in high resolution. Moreover, aliasing and flickering artifacts may occur because of the rasterization of the curves over a discrete multi-scale pixel grid. Jeschke et al. [2009a] present a faster solution, which initializes the color at each pixel with the color of the closest curve, and uses finite differences with variable step size to accelerate the convergence rate of Jacobi iterations. However, extensions that have been further proposed are not trivial to render with such an optimized approach.

These solvers are output sensitive, and thus possess important limitations: they become slow for high-resolution images, and they must be updated for most image operations. Moreover, a simple operation such as zooming-in implies many difficulties to take into account the constraints which are not visible. Jeschke et al. [2009b] adapt their 2D solver to deal with deformations due to texture mapping on 3D surfaces, using axis-aligned warping techniques. Their method is still limited regarding instancing (it must be recomputed for each instance) and complex deformations (it does not permit environment mapping). [START_REF] Pang | Fast rendering of diffusion curves with triangles[END_REF] take a different approach whereby DCs are converted into a triangle-mesh representation and colors are linearly interpolated over each triangle. Positive mean-value coordinates give the color at a vertex as a weighted sum of the colors of visible curves. To evaluate visibility, they employ an exhaustive sorting algorithm. An alternative consists in performing ray tracing on the GPU via a pixel shader [START_REF] Bowers | A ray tracing approach to diffusion curves[END_REF]]. The approach of [START_REF] Sun | Diffusion curve textures for resolution independent texture mapping[END_REF] avoids the explicit computation of visibility thanks to Green functions, which provide direct evaluation at individual points. Although the method is exact for closed curves, it only provides an approximation for open curves. Moreover, these three techniques are limited to harmonic interpolation, which lacks smoothness and control abilities.

Like [START_REF] Pang | Fast rendering of diffusion curves with triangles[END_REF], our solver makes use of an intermediate triangulation. However, it is based on higher-order patches and elegantly solves the Bi-Laplace equation using FEM while allowing for the largest set of artistic controls.

FEM for the Bi-Laplace equation

Even though recalling the huge literature on solving the biharmonic equation using FEM is out of the scope of this paper, it is worth noting that in the context of geometric modeling, similar fourth order PDEs often arise in the construction or preservation of high smoothness surfaces that are subject to zero and first order constraints. To achieve interactive performance, various discrete approximations that lack convergence have been employed (see for instance [START_REF] Botsch | On linear variational surface deformation methods[END_REF] for a comprehensive survey). To address the artifacts that may be produced by these approximations, [START_REF] Grinspun | Computing discrete shape operators on general meshes[END_REF] described a discrete shape operator that builds on non conforming elements. More recently, [START_REF] Jacobson | Mixed finite elements for variational surface modeling[END_REF] presented an appealing mixed FEM discretization which exhibits the same order of complexity than former approximations, but math- ematically more sound while permitting to seamlessly hanlde first order constraints. However, all these techniques have been developed specifically for linear meshes, while, as we show in this paper, color gradients require higher order elements.

∇uC\uC v l |vr v l |⊗ ⊗|⊗ ⊗ v g l |gr ⊗|gr ⊗|⊗ ⊗ g • | • • | ⊗ ⊗ | ⊗ ⊗ • 0|⊗ ⊗|⊗ ⊗ 0 (

Free-form vector gradients

The main goal of this paper is to produce smooth free-form gradients controlled by a sparse set of curve and point vector primitives. [START_REF] Finch | Freeform vector graphics with controlled thin-plate splines[END_REF] have shown that the solution to the Bi-Laplace equation ∆ 2 u = 0 provides for gradients that closely mimic smooth shading. It also offers advanced artistic controls in the form of user-specified value and gradient constraints, as illustrated in Figure 2. We thus target the same solution, considering u as an arbitrary scalar function (color channel, height, etc). In this section, we provide a taxonomy of constraints that unifies previous methods and further extends artistic controls. The way these constraints are handled in our solver is presented in Section 4.

Taxonomy of constraints.

For the sake of simplicity, in this paper we will assume curves are provided as Bézier splines, though any other representation could be used. At any point along a curve C, the value and/or the gradient magnitude in the normal direction may be constrained independently on each side of the curve. Denoting by uC the restriction of u to C, we summarize the space of possible constraints by the following diagrams:

= v l ,⊗ left v r ,⊗ = g l ,⊗ g r ,⊗ u c : right ∇ u c : left right
where l and r denote the 'left' and 'right' sides respectivelly. For instance, the right side of C may be assigned a specific scalar function vr or be left unspecified (we use the symbol ⊗ in such cases). Moreover, left and right values may be constrained to be equal (we then write uC = v). Similar choices are available for gradient constraints. The equality constraint is required whenever values or gradient magnitudes are left unspecified and one wants to enforce either C 0 or C 1 continuity. It also facilitates user interactions when the same values have to be set on both sides of a curve. As with previous work on DCs, scalar and gradient magnitude functions along curves are obtained through arc-length interpolation of an arbitrary number of constraint control points attached to curves at arbitrary locations. The case of points is slightly different: we treat them as isotropic constraints, and let users control not only their value and gradient magnitude, but also the gradient direction.

This small and intuitive set of constraints permits to achieve a large number of curve types. They are summarized in Figure 3(a), where we show all possible combinations of value and gradient constraints, each cell representing one type of curve. Such a taxonomy unifies and generalizes DCs and their extensions presented in earlier work, like barriers or slopes; these are directly identified by color in the table for completeness. In particular, observe that up to now, the use of gradient constraints has been limited to curves that have the same constraints on each of their sides (the 4 colored cells at bottom right). The offset curves of [START_REF] Finch | Freeform vector graphics with controlled thin-plate splines[END_REF] are not referenced since these are more a matter of user interface.

A benefit of our classification is that it lets emerge novel curve types that make more complex usages of gradient constraints (white cells), even though not all curve types make sense for all applications. For instance, choosing a specific gradient magnitude in color images makes little sense, appart from ∇uC = 0 that we use in some instances in the paper. In contrast, when creating height fields (terrains), a gradient constraint ∇uC = g is useful since it constrains the normal of the underlying surface. A selected subset of curve types is illustrated in Figure 3 Contours. When values or gradient magnitudes on one or both sides of a curve are left unspecified, additional constraints at the curve level may be enforced. This is already found in [START_REF] Finch | Freeform vector graphics with controlled thin-plate splines[END_REF], where they define contour curves as curves that have an unspecified, yet constant value along one side. In our taxonomy, such iso-values may be assigned independently to the right and/or left side of a curve, and take on arbitrary gradient constraints. Refering to the table in Figure 3(a), this applies to the three middle columns.Various configurations of iso-values are illustrated in Figure 4. Note that we could similarly enforce iso-gradients, which correspond to the three middle rows of the table.

Transmission. When values on both sides of a curve are left unspecified but enforced to be equal (⊗ column), any value in a family of solutions may be chosen to meet the equality constraint. We thus request an additional artistic parameter α ∈ [-1, 1] for such type of curves that permits to pick one specific solution. Intuitively, α controls the relative amount of diffusion across a curve: with α = 0, values are partially blocked on one side, but transmitted on the other; α = 0 boils down to a symmetric transmission; and α = ±1 performs a purely asymmetric transmission. As shown in Figure 5, such a transmission control is very useful since it permits to control the influence of other curves. Nil curves (gray cell in Figure 3(a)) are particularly useful when combined with transmission. As opposed to the method of [START_REF] Bezerra | Diffusion constraints for vector graphics[END_REF] that relies on global homogeneous color coordinates for this purpose, the α parameter provides an accurate local control over diffusion.

Vectorial solver

The main technical contribution of this paper is on the evaluation of the sparse representation defined in the previous section. This amounts to solving the Bi-Laplace equation ∆ 2 u = 0, subject to the constraints of Figure 3, 4 and 5. In contrast to previous approaches, we propose to output an intermediate denser vector representation enabling closed-form evaluation. We first present this representation, which consists in a set of triangular patches (Section 4.1). We then explain how we create it in two steps: the unknown function u is discretized into triangles covering the domain of resolution (Section 4.2), after which unknown patch parameters are solved using an appropriate Finite Element Method (FEM) (Section 4.3). This is depicted in Figure 6, where we also show the final result obtained after rendering (Section 5).

Intermediate representation

Central to our framework is the design of our intermediate vectorial representation. Quadtree-based meshes are fast to generate [START_REF] Frey | Fast adaptive quadtree mesh generation[END_REF]] but they would essentially lead to the same issues than directly solving in a pixel grid. Meshless representations [START_REF] Belytschko | Meshfree and Particle Methods[END_REF] are easy to manipulate and known to produce very smooth results, but they are not well suited to faithfully represent discontinuities and are relatively costly to evaluate. A triangular mesh with curved edges is therefore a natural choice since it is easy to make its vertices and edges match input constraints. More precisely we consider non-overlapping polynomial patches ui such that the function u is then represented by:

u(x) = u idx(x) (x) , (1)
where idx(x) returns the index of the unique triangular patch containing the evaluation point x. This choice significantly simplifies and speeds up our FEM solver, rasterization, and random evaluations. Using a higher-order representation such as subdivision surfaces (as in [START_REF] Liao | A subdivision-based representation for vector image editing[END_REF]) would greatly complexify these operations, hence limiting the benefits of an intermediate representation.

We use quadratic Lagrange polynomials: ui(x) = 6 k=1 v k B k (x) for our patches, where B k is the basis function associated to the k th node of the patch with associated value v k . Let Φ k , k = 1, 2, 3 be the linear basis functions associated with the triangle vertices (i.e., barycentric coordinates); then, using the notation of Figure 7a, we have B1 = Φ1(2Φ1 -1) at a vertex node, and B4 = 4Φ1Φ2 at an edge node. The other basis functions are obtained by cyclic permutations. As shown in Figure 8, quadratic patches outperform linear ones while remaining compact (6 coefficients) and fast to evaluate. Even though cubic polynomials might constitute an appealing alternative, they have not been considered because they would significantly increase the computational complexity of our FEM solver.

Representing singularities. Curve extremities or intersections might require to specify two different values at the same node. This creates a singularity that cannot be represented with standard triangular patches. Let's consider an extremity point p1 for instance, with value constraints v1→2 and v1→3, one for each abutting edge in the triangulation (see Figure 7b). Our key observation here is that at such an extremity, we expect the bi-harmonic interpolation to yield a smooth radial gradient from one edge to the other. To capture this behavior, we make the value v1 a function of the evaluation point x, performing a linear interpolation between v1→2 and v1→3 based on the relative angular position:

v1(x) = v1→2 + (v1→3 -v1→2) (x -p2) T (x -p1) (p3 -p2) T (x -p1) .

Constrained triangulation

The goal of the triangulation step is to produce a triangular partition of the domain that satisfies the following criteria:

• it must exactly match the input constraints,

• it must be dense enough to faithfully represent u,

• it must be fast enough to permit interactive drawing,

• it must be locally uniform to avoid numerical issues and provide high visual quality.

Note that the last three criteria may be in conflict: a high quality often requires a high density, which also demands high computational and memory costs. To achieve these goals, we implemented a fourstep triangulation procedure, illustrated in Figure 9. The first step consists in the construction of a planar straight lines graph (PSLG) by discretizing the input curves into a set of polylines. Then a constrained Delaunay algorithm yields an initial triangulation, which is refined by two additional steps: a refinement around singularities, and a geometric refinement that provides an overall high-quality (dense enough) triangulation. Our procedure may be seen as a variant of the work of [START_REF] Pang | Fast rendering of diffusion curves with triangles[END_REF], and fixes several shortcomings as detailed below.

Construction of the PSLG. The purpose of the first step is mostly to simplify and speedup the Delaunay triangulation and FEM solver by approximating input curves by polylines. In order to satisfy our first criterion, we construct the PSLG such that it includes all input constraint points as well as curve intersections.

The remaining curve segments are then subdivided until the maximal angle between the segment and tangent extremities is lower than a given threshold (15 • in our examples). It is important to discretize input curves in such an order, otherwise two vertices might be generated arbitrarily close to each other, thus leading to overrefinement in the final Delaunay refinement step.

Singularity refinement. The strongest variations often occur around curve singularities (tips and intersections) and isolated point constraints. To guarantee high quality around such a point p, we insert up to six vertices qi uniformly spread around p at a distance l such that the radial angle qipqi+1 between two successive vertices is smaller or equal to π/3. For tips or intersection points, the qi are inserted along the constrained edges in priority. The insertion distance l is adaptively determined from the Delaunay triangulation which provides meaningful information about the relative distances between neighboring constraints. In particular, we chose l as the third of the shortest triangle altitude passing through p (Figure 9c).

Geometry refinement.

Since quadratic polynomials cannot represent inflections, we must ensure that two constraints are not directly connected by an edge. In order to limit the number of generated triangles, we satisfy this requirement through Delaunay refinement [START_REF] Shewchuk | Mesh generation for domains with small angles[END_REF]] while ensuring a minimal triangle quality. Given a triangle metric, Delaunay refinement iteratively inserts a new vertex into the constrained Delaunay triangulation for the triangle having the worst score. To achieve our goals, triangles connecting two constraints are prioritized and then sorted according to the ratio between the circumradius and the shortest edge length, which is a standard quality metric implicitly bounding the minimal and maximal angles [START_REF] Shewchuk | Mesh generation for domains with small angles[END_REF]]. For our results, we stopped the refinement once the minimum triangle angle is greater than 20 • . We have found these a priori refinement heuristics preferable to more complex mechanisms that would violate our performance criterion. For instance, interleaving re-meshing and solving operations would likely provide a better control over density, but would be prohibitively expensive to compute. Regarding performance, it is important to remark that our triangulation procedure only depends on input geometry; hence it does not have to be recomputed when constraint values are modified while their positions are kept unchanged.

FEM solver

At this stage, we have discretized the function u into adjacent quadratic triangular patches ui. This section explains how to make use of Finite Element Methods (FEM) to accurately and efficiently compute nodal values vj such that u satisfies the Bi-Laplace equation.

As discussed in Section 2.4, similar problems have already been studied in the graphics community. However, all these techniques have been developed specifically for linear meshes, while, for our application, quadratic patches appear to be the best trade-off to convey visually pleasant results while keeping a coarse triangle mesh. Indeed, a general result in the FEM literature is that the convergence rate with respect to element order is higher than for increasing element resolution [START_REF] Zienkiewicz | The Finite Elements Method: Its Basis and Fundamentals[END_REF]]. As discussed in Section 4.1, we have confirmed that in our experiments: for color gradients quadratic patches offer an obvious higher visual quality than using four times more linear patches. The difference is amplified when editing or animating curves.

Moreover, our application exhibits special needs such as singularities and punctual gradient constraints, which are not straightforward to handle using the aforementioned discretizations. To address these difficulties, we propose to use a weak formulation of our differential equation as well as non-conforming elements as detailed below.

Weak formulations. FEM never directly solves for the strong formulation of a differential equation. The equation is rather projected onto an appropriate set of test functions tj in order to reduce requirements on continuity and degree of the elements, as well as to improve the stability of the numerical scheme [START_REF] Zienkiewicz | The Finite Elements Method: Its Basis and Fundamentals[END_REF]]. In our case, this leads to the following problem [START_REF] Lascaux | Some nonconforming finite elements for the plate bending problem[END_REF]:

∀j, Ω ∆u • ∆tj -σ ∂ 2 u ∂x 2 ∂ 2 tj ∂y 2 -2 ∂ 2 u ∂x∂y ∂ 2 tj ∂x∂y + ∂ 2 u ∂y 2 ∂ 2 tj ∂x 2 = 0 , ( 2
)
where σ is Poisson's coefficient, which controls the influence of the regularization terms. In practice, σ has to be chosen in the range [ 1 2 , 1], in which case its value does not influence the result. Choice of elements. Ideally, we would like to directly use our current discretization of u in quadratic triangular patches to solve Equation 2. Unfortunately, the numerical scheme requires C 1 continuity across such conforming elements to reach convergence. This typically implies the use of quintic polynomial elements, which are prohibitively expensive to compute.

w j w j - w j w k v = •/g = • v = •/g = •|• v = •| • /g = • v = •| • /g = •|•
Non-conforming elements permits to reach convergence without requiring C 1 continuity across triangle edges. For instance, the Morley element [START_REF] Morley | On the constant moment plate bending element[END_REF] shown in Figure 7c is the lowest-degree element that can solve Equation 2. It is a quadratic polynomial basis function with six degrees of freedom: the three nodal values at the vertex positions, and the three normal derivatives at the edge middles. Note that a variant of Morley's elements have been successfully used to obtain high quality and stable solutions for geometric modeling purposes [START_REF] Grinspun | Computing discrete shape operators on general meshes[END_REF]]. In practice, Morley's elements are not even C 0 continuous across their boundaries. C 0 continuity could be enforced by computing the nodal edge values of our Lagrange patches from the average of the edge middle values of the two adjacent Morley elements. Unfortunately, as shown Figure 8c this naive strategy leads to very unpleasant results.

In order to be able to output quadratic patches with at least C 0 continuity, we thus propose to use the third-order non-conforming Fraeijs de Veubeke's element [1974] (FV) which is illustrated in Figure 7d. The FV element can be seen as a higher degree variant of Morley's element. It is controlled by six nodal values at vertices and mid-points of the edges, and the mean value of the normal derivative along each edge (which we call mean gradient nodes). This element is perfectly suited to our problem. Indeed, since its nodal parameters form a super-set of parameters for our quadratic patches, the conversion from FV elements to our Lagrange polynomials is straightforward: we simply dismiss the mean gradient nodes. Moreover, it provides a control over the gradient magnitude in the direction normal to an edge, which is necessary for handling gradient constraints.

Under these settings, we are now seeking for an approximation of u of the form j wjFj where the Fj are the basis functions of the FV elements given in appendix, and the wj are the sought-for nodal values (with wj = vj for j ≤ m, m being the total number of nodes in our intermediate representation). As it is often done in FEM, we take for the test functions tj the same set of basis functions Fj. Then Equation 2 becomes a sparse linear problem of the form Aw = 0 where the matrix A is called the stiffness matrix whose coefficients Ai,j are given in the appendix.

In Section 4.4, we give some implementation details; in particular we explain how simple edits of constraint values only require partial updates of the system. In the rest of this section we show how to handle our specific constraints.

Handling constraints. Thanks to the use of FV elements, both value and gradient constraints are easily satisfied (because all constraint points and curves are vertices and edges of the triangulation).

As illustrated in Figure 10, the non-equality of values (resp. gradient) across a curve is handled by doubling the value (resp. mean gradient) nodes along the curves. In table 3 this concerns curves of the form •|•. A node attached to a constraint providing a specific value directly takes its value from the constraints, and it is removed from the list of unknowns (this corresponds to columns v l |vr, v l |⊗, and v). Likewise, gradient values specified on a curve are integrated over the respective edge of the element to fix the respective mean gradient node (this corresponds to rows g l |gr, ⊗|gr, and g). Isocontours are trivially handled by having a single nodal value, i.e., a single unknown, for all the nodes lying on the same iso-contour.

Taking into account the gradient specified at an isolated constrained point is slightly more difficult because the FV element only offers a control over the gradient magnitude along an edge. Let us consider the FV element µ(x) = 9 k=1 F k (x)w k as in Figure 7d, assuming the vertex p1 is an isolated constraint point with gradient g. Let e2 = p2 -p1 and e3 = p3 -p1 be the two incident edges, and µ2(t), µ3(t) the value of the element over these two edges. Constraining ∇µ(p1) = g, is equivalent to constraining [µ 2 (0) µ 3 (0)] T = ḡ, where ḡ = [e2 e3] -1 g is the projection of g into the frame spanned by e2,e3. Looking at the FV's basis functions (Equation 3), one can remark that the mean gradient node of a given edge does not influence the value along its edge. Therefore the derivative constraint in each direction can be enforced independently by constraining the mean gradient nodes w7 and w9 respectively. For instance, this yields for w7:

w7 = ḡ0 -k =7,9 (F 2 k ) (0)w k (F 2 7 ) (0)
, where F 2 k = F k (p1 + te2) are the FV basis functions restricted to the edge e2. These constraints are introduced in equation 4 when assembling the stiffness matrix, thus effectively removing two degrees of freedom.

Handling transmission.

Transmission permits values coming from one side of a curve to influence the other side, and may be made asymmetric by partly blocking the diffusion in the opposite direction. This is particularly useful to control highlights for instance, as seen in Figure 11. Such a curve has no value constraint except that values must be equal on each side (⊗ column). Observe that the unknown wj of a given node belonging to such a curve has an influence on the nodal values of the incident elements. Given two weights β l and βr, the key idea is to make the node j emit β l wj to the left nodes, and βrwj to the right nodes. When building the stiffness matrix, this amounts to inserting β l Ai,j instead of Ai,j if the node i is on the left, βrAi,j if the node i is on the right, and directly Ai,j otherwise (i is on the curve). Note that this leads to an asymmetric stiffness matrix. The weights β l and βr are computed from the transmission factor α. When α = 0, we set both β l = βr = 1, leading to the default behavior. With α > 0, βr = 1 and β l = 1 -α, thus elements on the right will behave normally, but the ones on the left will be less influenced by wj and therefore less influenced by the right side. Ultimately, when α = 1, β l = 0, elements on the left do not depend on values on the right side anymore, which effectively amounts to block diffusion from right to left. Conversely, when α < 0 we take β l = 1 and βr = 1 + α. Handling singularities. As shown in section 4.1, curve extremities and intersections require special patches to appropriately represent the singularity that appears when two different values are assigned to a same node. While our novel singular patch works well for evaluation purposes, it cannot be used for solving purpose since it is not differentiable. We workaround this issue by approximating a singular element with a pair of standard FV overlapping elements that are connected to two different nodes at the singularity. As shown in Figure 12, this simple approach works remarkably well in practice.

Implementation details

Our current implementation entirely runs on the CPU, and does not exploit multi-threading yet. It makes use of CGAL for the Delaunay triangulation and refinement, and Eigen's direct solver for solving our sparse problems. In contrast to iterative methods, a direct solver permits to factorize the left hand side matrix once, to then very efficiently solve for multiple right sides. For color images, this permits to efficiently solve for the red, green, blue, and alpha components at once. Moreover, since the modification of a constraint value only changes the right hand side of the system, the factorization does not have to be recomputed during constraint value editing. This implementation permits the setting of a rich number of constraints (like gradient magnitudes, contour curves or asymmetric transmission) directly in the linear system of equation. We have tested our method on one single core of a Intel Core i7-3610QM, and we report our performances in Table 1. It is difficult to compare our performances with those of the pixel-based method of [START_REF] Finch | Freeform vector graphics with controlled thin-plate splines[END_REF]. Indeed, the complexity of their solver grows with image resolution, while ours grows with the number of vector primitives. We believe that output insensitive solvers are preferable since the entire goal of vector drawing is to keep a small number of sparse primitives. Moreover, let us note that for images composed of a number of primitives similar to their results, we obtain similar performances irrespective of image resolution. On top of that, we have used a single core implementation, while they explicitly mention using a parallelized implementation on a 2.67GHz dual quad-core Intel Xeon. Our approach thus proves to be efficient for the interactive editing and drawing of free-form vector gradients. However, its main benefit is to output a vectorial representation that makes subsequent rendering tractable, as explained in the next section.

Rendering

Our solver needs to be fully updated everytime the position of a control point is modified, and partially when constraints are changed. Such modifications are part of the drawing/editing process. This is to contrast with rendering, which essentially consists in evaluating our intermediate representation at specific locations. Since our solution is given in closed-form, it tremendously simplifies and improves applications that make use of vector gradients, as detailed in the following.

Color images. With our approach, color images are best displayed by directly rasterizing the quadratic patches. Even though we assumed linear edges when solving, at rendering time edges tied to an input curve must be properly refined according to the zoom and deformation levels. In our implementation, such triangles are subdivided in software, and then rasterized using OpenGL and pixel shaders to evaluate the quadradic patches. Hardware tessellation could be used for this purpose too.

Our approach permits to naturally handles a broad range of deformations at rendering time: translations, rotations and scaling, but also perspective deformation as shown in Figure 1. One simply has to deform the intermediate representation itself. This is particularly useful for printing since arbitrary high image resolutions may be output (Figure 13a); for copy-pasting since copies are treated as (possibly deformed) instances with no increase in complexity (Figure 13b); and for storage and dissemination. In contrast, pixelbased solvers must impose additional constraints at image borders when zooming in; they must be fully updated for most deformations; and their complexity quickly increases with higher image resolutions and with additional deformed instances.

2.5D images.

Although not all of the gradient constraints from our taxonomy make sense for color images, they reveal their usefulness when working with 2.5D data such as height fields or normal maps. This is shown in Figure 14 with the construction of a height field. Recall our patches are only C 0 continuous across them. Even though this has been proven to be enough for the display of color gradients, the shading applied to a height-field might reveal the lack of C 1 continuity. Higher quality could be obtained at rendering time by, e.g., applying quadratic normal interpolation [START_REF] Vlachos | Curved pn triangles[END_REF]. As shown in Figure 15a, normal images may also be directly created with our system by propagating surface normal coordinates. In this case, the user does not set coordinate constraints on curves though: we enforce the projection of a 3D normal constraint to be aligned with the 2D normal of the curve as in [Winnemöller et al. 2009]. This leaves a useful degree of freedom, which corresponds to the normal coordinate in the direction perpendicular to the image plane. Moreover, assuming quadratic surfaces, the normal field is expected to vary linearly in the image. We thus use a simplified version of our solver that performs harmonic interpolation (see appendix). Our approach has a number of benefits compared to previous methods for the creation of height fields [Winnemöller et al. 2009;[START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF]]. First, it provides a richer set of artistic controls thanks to our novel taxonomy. Second, it only requires a single optimization, while others require integration in post-process [Winnemöller et al. 2009] or multiple diffusion steps [START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF]]. Third, it outputs results directly in a structure amenable to rendering, while others need either a specific renderer or a mesh conversion in post-process. Indeed, our intermediate mesh only has to be lifted based on computed height values, and then be loaded in a rendering engine.

3D mapping. We have also adapted our method to the mapping of vector gradients onto parametrized surface meshes in 3D. One solution consists in intersecting the triangles of our mesh, with the triangles that define the surface in texture space, hence refining the 3D mesh to precisely meet gradient requirements. This is illustrated for color textures in Figure 16, and could be easily adapted to work with displacement maps. In this later case, the output mesh would naturally preserve the discontinuities (cliffs) present in the vector displacement map. Another solution consists in accessing randomly our intermediate representation at required texture coordinates. For instance, the realtime GPU based techniques of [START_REF] Parilov | Real-time rendering of textures with feature curves[END_REF], or [START_REF] Nehab | Random-access rendering of general vector graphics[END_REF] can already handle quadratic curves, and would require little effort to be adapted to our purpose. Again, our method provides a lot more flexibility compared to previous work. Instancing is a very common practice in texture mapping, and our solution does not require any specialized solver like [Jeschke et al. 2009b]. Thanks to random access evaluation, it works for arbitrary parameterizations, even environment maps as shown in Figure 15. To the best of our knowledge, no previous method has shown such uses of free-form vector gradients.

Discussion

Summary. The vectorial solver presented in this paper combines for the first time the expressivity of sparse curve and point vector drawing primitives, with the rendering simplicity of vector representations. This is permitted by an intermediate triangular structure that is updated interactively during drawing, but most importantly that is trivially evaluated during rendering, even at random locations. Our method unifies and generalizes the types of controls achievable by previous techniques, and introduces a novel trans- mission parameter that provides a much-needed local control over propagation. The use of FEM with high-order elements permits an accurate and fast solution to the Bi-Laplace equation, and directly deals with gradient constraints. In contrast, other solutions are either limited to the Laplace equation [START_REF] Orzan | Diffusion curves: A vector representation for smooth-shaded images[END_REF][START_REF] Bezerra | Diffusion constraints for vector graphics[END_REF]Jeschke et al. 2009a], simple gradient constraints [START_REF] Finch | Freeform vector graphics with controlled thin-plate splines[END_REF], or require multiple passes with indirect gradient control [START_REF] Bezerra | Diffusion constraints for vector graphics[END_REF][START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF].

Limitations. Triangulation updates might lead to flickering artifacts when drawing primitives are edited; although this is a major problem with linear elements, our quadratic elements reduce this issue to a point where they are hardly noticeable, as is best observed in the accompanying video. Even though quadratic elements are more costly for optimization, they demand significantly simpler triangulation to reach a reasonable quality and are thus comparable in terms of performance to linear elements. Our optimization procedure could be improved in a number of ways: employing a parallelized implementation, identifying closed curves to reduce dependencies, or using a cheaper triangulation than Delaunay (care should be taken with degenerated triangles though). Some limitations of our system are not due to the solver itself, but rather to the choice of the Bi-Laplace equation. First, as identified by [START_REF] Finch | Freeform vector graphics with controlled thin-plate splines[END_REF], when only colinear constraints are used, the solution is not determined and we must revert in this case to the Laplace equation. But more importantly, the Bi-Laplace equation may create new value extrema outside from curve or point constraints, whereas the Laplace equation is guaranteed to avoid such issues. This is especially noticeable when using gradient constraints at singularities. This isssue has been recently addressed through a non linear optimization guided from an harmonic solution [START_REF] Jacobson | Smooth shape-aware functions with controlled extrema[END_REF]]. An alternative solution might be to blend between Laplace and Bi-Laplace equations in problematic regions. We believe we could adapt our solver to this purpose in future work since it is already able to solve Laplace equation.

Future work.

In this paper we have focused on free-form vector drawing, but our approach may also be used for image vectorization. The goal would then be three-fold: identify vector primitives and their type in an image; generate a triangulation with a density adapted to remaining variations; find color and gradient constraints that best approximate the input image inside each triangle. Contrary to previous techniques (e.g., [START_REF] Liao | A subdivision-based representation for vector image editing[END_REF]) this approach would attach constraints to a few vertices, making subsequent manipulations a lot more direct since the interior triangulation could be safely modified during editing. We have also restricted our solver to work in a flat 2D domain. It could be adapted to more complex domains in future work, for the drawing of free-form gradients directly on surfaces, or for creating volumetric gradients in 3D.

Figure 2 :

 2 Figure 2: Simple curves. The top and bottom curves diffuse color gradients on both images. The middle curve either acts as a diffuser with different colors on each side (left), or as a barrier (right).

  a) Taxonomy of constraints (b) Variations on color gradients

Figure 3 :

 3 Figure 3: Types of curves. (a) By combining value constraints (columns) and gradient constraints (rows), we define a rich set of curve types. Their specificities are best understood by inspecting them row by row, observing how different types of value constraints yield new configurations (points & dotted lines correspond to value & gradient constraints). Our taxonomy generalizes previous work, as indicated by cell colors that correspond to existing types of DCs: diffusion, barrier, tear, crease, crease-value, value, value-slope, slope. The nil curve has no effect since it sets no constraints. (b) The four last rows of the table are illustrated with color gradients. The top and bottom curves of Figure 2 are left unchanged, but the type of middle curve is modified. We overlay iso-values of primary colors for visualisation purposes.

Figure 4 :

 4 Figure 4: Contour curves may be used with different types of value constraints, as shown here for the middle curve. Images from left to right correspond to ⊗ , ⊗|vr, and ⊗|⊗ columns.

Figure 5 :

 5 Figure 5: Transmission curves permit to locally control color diffusion, here using the nil curve in the middle. From left to right, we have used α = -1, α = 0, and α = 1.

Figure 6 :

 6 Figure 6: Overview of our solver. Input vector primitives (a) are first discretized to produce our intermediate representation (b). It is used as input to our FEM solver that outputs a set of elements (c), which are easily rendered to obtain the final image (d).

Figure 7 :

 7 Triangular patches and elements. (a) A Lagrange patch with its six nodes. (b) Our special quadratic patch for singularities. (b) Morley's element. (d) Fraiejs de Veubeke's element.

Figure 8 :

 8 Figure 8: Patch order. We provide a visual comparison between 2400 linear (a), and 1200 quadratic (b) patches. Both images requires about 5000 nodes, and have been computed in 60ms each. On the bottom, grey levels are mapped to RGB values to better reveal artifacts. Converting Morley's element to quadratic patches (c) leads to even stronger artifacts.

Figure 9 :

 9 Constrained triangulation. (a) Curve distretization yields a PSLG in green. The other figures show a closeup view of the triangulation process: (b) Delaunay triangulation; (c) Explicit refinement around critical points; (d) Geometric refinement.

Figure 10 :

 10 Figure 10: Constraints along a curve are obtained by sharing or duplicating value or gradient nodes.

Figure 11 :

 11 Figure 11: Transmission curve (dashed) is used to create a complex highlight and blended with global gradients (plain curves).

Figure 12 :

 12 Figure 12: Singularities are handled by duplicating elements (left). Results obtained with a coarse or very dense triangulation are indistinguishable (right). The highlighted triangles correspond to the singular ones with in bold the one of the left figure.

  Fig 1a Fig 8 Fig 11 Fig 13b

Figure 13 :

 13 Color results. Our approach is independent of image resolution (a), and is especially efficient for instancing (b).

Figure 14 :

 14 Figure 14: Height field construction. A few curves (top right) suffice to create a terrain. Here, green curves impose value constraints, while red curves set gradient constraints only.

  Texture mapping. The texture mapping shown in Figure 1 is obtained by: (a) positioning a vector color image in texture space; and (b) intersecting surface triangles with triangles coming from our representation (shown before and after intersection).

Table 1 :

 1 Performances (in ms) for the different steps of our solver on a few examples. Images with asymmetric transmission yield nonsymmetric sparse problems which are more involved to solve.

		(symm.)	(non-symm.)	(non-symm.)	(symm.)
	PSLG	41	2	1	8
	Refinement	38	12	3	9
	Assembly	65	16	7	14
	Solve	58	35	14	13
	Total	202	65	25	44

This work has been supported by the ANR SeARCH (ANR-09-CORD-019) and ANR DRAO projects. We would like to thank Noémie-Fleur Sandillon-Rezer for designing Figures 1 and13a, and Laurence Boissieux for the model of Figure 6.

Appendix

Fraeijs de Veubeke's element. For the sake of reproducibility, the closed form formulas of the basis function for the FV element are given below:

. Here, Φi corresponds to the linear basis elements (barycentric coordinates), ∆ is the area of the triangle, ei is the edge opposite to the i th vertex. The quantities d2, d3, and the other basis are obtained by circular permutation. The coefficient Ai,j of the stiffness matrix are of the form:

Harmonic interpolation Solving the Laplace equation ∆u = 0 with our approach is considerably simpler. Its associated weak formulation ∀j, Ω ∇u • ∇tj = 0 , requires only C 0 continuity across elements. Therefore, it can be directly solved using our Lagrange patches as the elements. On the other hand, gradient constraints cannot be imposed in this case.