Michael Döhler
email: michael.doehler@inria.fr

Qinghua Zhang
email: qinghua.zhang@inria.fr

Laurent Mevel
email: laurent.mevel@inria.fr

Vibration Monitoring by Eigenstructure Change Detection Based on Perturbation Analysis

Keywords: Vibration monitoring, eigenstructure, statistical tests, fault detection

Vibration monitoring, notably in the fields of civil, mechanical and aeronautical engineering, aims at detecting damages at an early stage, in general by using output-only vibration measurements under ambient excitation. In this paper, a new method is proposed for the detection of small changes in the eigenstructure of such systems. The main idea is to transform the multiplicative eigenstructure change detection problem to an additive one, by means of perturbation analysis based on the assumption of small eigenstructure changes. Another transformation then further simplifies the detection problem into the framework of a linear regression subject to additive white Gaussian noises, leading to a numerically efficient solution of the considered problem. Compared to existing methods, it has the advantages of focusing on chosen system parameters and efficiently addressing random uncertainties. A numerical example of a simulated mechanical structure and a lab experiment on a beam, each with the detection of different damages, are reported.

INTRODUCTION

The detection of damages based on measured vibration data is a fundamental task for structural health monitoring to allow an automated damage diagnosis, e.g., for mechanical, civil or aeronautical structures [START_REF] Farrar | An introduction to structural health monitoring[END_REF]. If damages are detected early enough, they can be accommodated by maintaining the structures in appropriate operational states. Damage affects the dynamic properties of a structure, inducing changes in the modal parameters (natural frequencies, damping ratios, mode shapes), or in their equivalent eigenstructure representation (eigenvalues and observed eigenvectors). Damage detection thus amounts to detecting small changes in the eigenstructure of the monitored structure. A particular difficulty for structural health monitoring is caused by the absence of known system inputs, since the structural excitation is usually only ambient, leading to an outputonly monitoring problem.

There exist various methods for mechanical structure damage detection [START_REF] Carden | Vibration based condition monitoring: a review[END_REF], either modelbased or data-based. Model-based detection methods typically fall into two categories, respectively based on repeated system identification and on direct model-data matching.

Within the methods based on repeated system identification, a common strategy is to repeatedly estimate current modal parameters by means of system identification, and to compare the result to some reference modal parameters [START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: From automated OMA to damage detection[END_REF]. Tracking the evolution of param-This work has been supported by the ITEA2 MODRIO project.

eters over time, recursive identification schemes may be used, e.g., based on subspace methods [START_REF] Goethals | Recursive output-only subspace identification for inflight flutter monitoring[END_REF], on maximum likelihood estimation using a Kalman filter coupled with a tangent filter [START_REF] Campillo | Recursive maximum likelihood estimation for structural health monitoring: tangent filter implementations[END_REF] or on interacting Kalman filter and particle filter [START_REF] Zghal | Modal parameter estimation using interacting Kalman filter[END_REF].

With the methods based on direct model-data matching, current measurement data are directly confronted to a reference model, without resorting to repeated system identification. For instance, such methods include nonparametric change detection based on novelty detection [START_REF] Worden | Damage detection using outlier analysis[END_REF][START_REF] Yan | Structural damage diagnosis by Kalman model based on stochastic subspace identification[END_REF] or whiteness tests on Kalman filter innovations [START_REF] Bernal | Kalman filter damage detection in the presence of changing process and measurement noise[END_REF]. Another method within this category, the local asymptotic approach to change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF], has the ability of focusing the detection on some chosen system parameters, in particular on modal parameters of mechanical structures. Associated to efficient hypothesis testing tools, this method has led to successful applications in the field of vibration monitoring [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF], [START_REF] Zouari | Subspace-based damping monitoring[END_REF], [START_REF] Jhinaoui | Subspace instability monitoring for linear periodically time-varying systems[END_REF], [START_REF] Döhler | Subspace-based fault detection robust to changes in the noise covariances[END_REF], Döhler et al. [2014a,b]).

The purpose of the present paper is to propose another method based on direct model-data matching, which can also focus the detection on chosen system parameters, and has the advantage of more efficiently addressing random uncertainties when associated to hypothesis testing tools. In this method, the normalization of the involved statistical tests is fully based on the covariance matrices of the state and output noises estimated at the step of the nominal system identification, whereas in the local asymptotic approach to change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF], the statistical tests are normalized by the covariance matrix of the sum of some designed residual vector, whose estimation requires large data samples.

The eigenstructure change detection problem is multiplicative in the sense that the eigenstructure parameters appear in coefficients of the unknown state vector of the state space model characterizing mechanical structures. The main idea of the present paper is to transform this problem to an additive change detection problem, by applying the perturbation analysis based on the assumption of small eigenstructure changes. After this transformation, another step based on the results of [START_REF] Zhang | Statistical detection and isolation of additive faults in linear time-varying systems[END_REF]) then further transforms the problem into the form of a linear regression with additive Gaussian noises. The covariance matrix of these Gaussian noises is easily determined from the properties of the nominal state space model of the monitored system, typically obtained by means of system identification.

This paper is organized as follows. In Section 2, the structural modelization is recalled and the fault detection problem is stated. In Section 3, the perturbation analysis is carried out to transform this problem into an additive one, and in Section 4 the respective hypothesis test is derived. An algorithmic summary of the developed method is given in Section 5. Finally, applications for vibration-based damage detection are shown in Section 6 in simulations and on a lab structure.

PROBLEM STATEMENT

The behavior of linear time-invariant mechanical structures subject to unknown ambient excitation can be described by the differential equation M Ẍ (t) + C Ẋ (t) + KX (t) = υ(t)

(1) where t denotes continuous time; M, C, K ∈ R m×m are mass, damping, and stiffness matrices, respectively; the state vector X (t) ∈ R m is the displacement vector of the m degrees of freedom of the structure; and υ(t) is the external unmeasured force (noise).

Observed at r sensor positions (e.g. displacement, velocity or acceleration sensors) at discrete time instants t = kτ (with sampling rate 1/τ ), system (1) can also be described by a discrete-time state space system model [START_REF] Juang | Applied system identification[END_REF]

z k+1 = F z k + ω k y k = Hz k + v k (2)
where the state vector z k = X (kτ ) T Ẋ (kτ ) T T ∈ R n with n = 2m, the measured output vector y k ∈ R r , the system matrices

F = exp 0 I -M -1 K -M -1 C τ ∈ R n×n H = L d -L a M -1 K L v -L a M -1 C ∈ R r×n , with selection matrices L d , L v , L c ∈ {0, 1}
r×m indicating the positions of displacement, velocity or acceleration sensors. The state noise ω k and output noise v k are unmeasured and assumed to be Gaussian, zero-mean, white.

Damage leads to changes in the structural properties of system (1), e.g., in mass parameters of elements in M, or in parameters corresponding to element stiffness such as Young's modulus in K. Hence they provoke changes in the eigenstructure of system (1), and consequently of system (2), which shall be monitored for changes. The eigenstructure of ( 2) is defined as the collection of eigenvalues and observed eigenvectors (λ i , ϕ i ) with

F φ i = λ i φ i , ϕ i = Hφ i , i = 1, 2, . . . , n (3)
and constitutes a canonical parameterization invariant to linear state transformations.

Assume that the eigenstructure of the considered system contains only complex modes. This is the typical case for structural health monitoring applications. Then, the eigenvalues λ i of the real matrix F ∈ R n×n (n = 2m) consist of m conjugate complex pairs. Let the vector

λ [λ 1 , λ 2 , . . . , λ m ] T ∈ C m contain m of the n eigenvalues λ i ,
one out of each of the m conjugate complex pairs, and ϕ [ϕ 1 , ϕ 2 , . . . , ϕ m ] ∈ C r×m be composed of the corresponding observed eigenvectors. The complex eigenstructure (λ i , ϕ i ) is then represented by the equivalent real parameter vector θ ∈ R n+nr defined as

θ    Re(λ) Im(λ) Re(vec(ϕ)) Im(vec(ϕ))    (4)
where Re and Im denote respectively the real part and the imaginary part of a complex value.

Assume that the matrix F in (3) is diagonalizable. Let the real matrices A(θ) ∈ R n×n , C(θ) ∈ R r×n be defined as

A(θ) = Re(diag(λ)) Im(diag(λ)) -Im(diag(λ)) Re(diag(λ)) (5a) C(θ) = [Re(ϕ) Im(ϕ)]
(5b) then the state space system (2) is equivalent to

x k+1 = A(θ)x k + w k (6a) y k = C(θ)x k + v k (6b
) in the sense that there exists an invertible matrix T ∈ R n×n such that x k = T z k and w k = T ω k . This eigencanonical state space model will be used in this paper for detecting changes in the eigenstructure.

In practice, the mass, damping, stiffness matrices M, C, K ∈ R m×m are usually unknown. To characterize the considered system, a state space model of order n = 2m can be estimated from available output data y k by methods of system identification [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation[END_REF], [START_REF] Ljung | System Identification -Theory for the User[END_REF]). The computation of the eigenvalues and eigenvectors of the estimated system matrix then leads to the nominal parameter vector, say

θ 0     Re(λ 0 ) Im(λ 0 ) Re(vec(ϕ 0 )) Im(vec(ϕ 0 ))     , (7)
characterizing the nominal state space model in the eigencanonical form as defined by ( 5) and ( 6). The covariance matrices of the state and output noises,

Q S S T R = E w k v k w T k v T k , (8)
are also derived from the result of system identification and from the transformation to the eigen-canonical form.

Assume that the currently monitored system is subject to a small change of the parameter vector, so that the value of θ becomes θ = θ 0 + εθ 1 (9)

where the unknown vector θ 1 ∈ R n is expressed as

θ 1 =     Re(λ 1 ) Im(λ 1 ) Re(vec(ϕ 1 )) Im(vec(ϕ 1 ))    
and ε > 0 is a small constant reflecting the fact that the parameter increment θ εθ 1 (10) is small. Under the small change assumption, the parameter vector θ deviated from the nominal value θ 0 remains approximately an eigenstructure parameter vector.

The remaining part of this paper is for the purpose of developing a method for detecting small changes θ = εθ 1 possibly affecting the parameter vector θ.

PERTURBATION ANALYSIS

The main difficulty for detecting changes in the parameter vector θ resides in the fact that the matrices A(θ), C(θ), possibly subject to unknown changes of θ, appear in products with the unknown state vector x k in the state space model ( 6). The purpose of this section is to transform this multiplicative change detection problem into an additive form by means of perturbation analysis.

As the matrices A(θ), C(θ) defined in (5) are linearly parametrized by the vector θ as composed in (4), the decomposition of θ in (9) leads to

A(θ) = A(θ 0 ) + εA(θ 1 ) C(θ) = C(θ 0 ) + εC(θ 1 ).
For more compact notations, let us define

A 0 A(θ 0 ), A 1 A(θ 1 ) C 0 C(θ 0 ), C 1 C(θ 1 ) then A(θ) = A 0 + εA 1 (13a) C(θ) = C 0 + εC 1 . ( 13b
)
Based on the nominal system model characterized by θ 0 , the Kalman filter can compute the one-step-ahead state prediction, denoted by the short notation

x 0 k x k|k-1 (θ 0 ).
Le K be the steady state Kalman gain, then x 0 k is recursively computed as

x 0 k+1 = A 0 x 0 k + A 0 K(y k -C 0 x 0 k ) (14a) x 0 0 = E(x 0 ) (14b
) where E(x 0 ) is the mean value of x 0 based on prior knowledge. Because the monitored system is possibly subject to changes, the predicted state x 0 k based on θ 0 may be biased. Assume that

x k = x 0 k + εx 1 k (15)
where x 1 k ∈ R n is an unknown vector such that εx 1 k corresponds to the small bias caused by the small parameter increment θ = εθ 1 .

In the state equation (6a), substitute A(θ) with (13a), then

x k+1 = A 0 x k + εA 1 x k + w k
In the term εA 1 x k , replace x k with (15), then

x k+1 = A 0 x k + εA 1 x 0 k + ε 2 A 1 x 1 k + w k Omit the term involving ε 2 (because ε is small), it then yields x k+1 ≈ A 0 x k + εA 1 x 0 k + w k . ( 16
)
Remind that A 1 = A(θ 1 ) is linearly parametrized by θ 1 , with A(•) as defined in (5). To make explicit the linear dependance on θ 1 of the term εA

1 x 0 k , divide x 0 k into two sub-vectors x0 k ∈ R m and x 0 k ∈ R m such that x 0 k = x0 k x 0 k , then εA 1 x 0 k = εA(θ 1 )x 0 k = ε Ψ k θ 1 with Ψ k diag(x 0 k ) diag(x 0 k ) 0 m×nr diag(x 0 k ) -diag(x 0 k ) 0 m×nr , (17)
following from (5a). The (approximate) equation ( 16) is then rewritten as

x k+1 ≈ A 0 x k + Ψ k θ + w k (18)
with the parameter increment θ as defined in (10).

In equation ( 18), the matrix Ψ k is composed of sub-vectors of the estimated state x 0 k , which is computed with the Kalman filter ( 14), hence θ is the only unknown in the term Ψ k θ. Therefore, the parameter increment θ appears additively in the state equation ( 18), in contrast to its multiplicative form in (6a).

The term C(θ)x k in (6b) can be addressed similarly, thus

y k = (C 0 + εC 1 )x k + v k = C 0 x k + εC 1 x k + v k = C 0 x k + εC 1 (x 0 k + εx 1 k ) + v k ≈ C 0 x k + εC 1 x 0 k + v k
where again the term involving ε 2 was omitted.

In the last equation, since C(θ) is linearly parametrized as in (5b),

εC 1 x 0 k = εC(θ 1 )x 0 k = ε Φ k θ 1 = Φ k θ with Φ k 0 r×n (x 0 k ) T ⊗ I r . ( 19
)
and θ as defined in (10).

To summarize, the initially considered multiplicative change detection problem formulated with the state space model ( 6) is now transformed to the detection of additive changes in the new state space model

x k+1 ≈ A 0 x k + Ψ k θ + w k (20a) y k ≈ C 0 x k + Φ k θ + v k (20b)
where y k is the currently measured output, A 0 , C 0 are known from the nominal model (typically obtained by means of system identification), Ψ k , Φ k are composed of x 0 k which is estimated from the currently measured output y k through the Kalman filter (14).

ADDITIVE CHANGE DETECTION

Additive change detection problems in the form of ( 20) have been studied in [START_REF] Zhang | Statistical detection and isolation of additive faults in linear time-varying systems[END_REF]), by transforming the dynamic system model ( 20) into an equivalent linear regression model through a particular Kalman filtering.

Apply the Kalman filter to the state space system (20) while assuming θ = 0, it yields the one-step-ahead state prediction x k|k-1 (θ 0 ). Denote the innovation sequence (the prediction error) of this Kalman filter as

ζ k y k -C 0 x k|k-1 . ( 21
)
If θ = 0, (actually no eigenstructure parameter change), then it is well known that the innovation sequence ζ k is a centered Gaussian white noise.

However, if θ = 0, the innovation sequence ζ k is biased, because it is computed with the Kalman filter assuming θ = 0. In this case, according to the Proposition 2 of [START_REF] Zhang | Statistical detection and isolation of additive faults in linear time-varying systems[END_REF]), the innovation sequence ζ k satisfies (the "≈" would be replaced by "=" if the equalities in (20) were accurate)

ζ k ≈ (C 0 Γ k + Φ k ) θ + e k (22) where Γ k is recursively computed as Γ k+1 = A 0 (I n×n -KC 0 )Γ k + Ψ k -A 0 KΦ k (23a) Γ 0 = 0, ( 23b
) and e k is a white Gaussian noise of zero mean. The covariance matrix of e k is Σ = C 0 P (C 0 ) T + R (24) where P is the solution of the algebraic Riccati equation associated to the Kalman filter and R is the covariance matrix of the output noise v k .

It turns out that this Kalman filter delivering the state prediction x k|k-1 and the associated innovation sequence ζ k in ( 21) is exactly the same as the one expressed in ( 14), as both are based on the same nominal system matrices A 0 , C 0 , the same noise covariance matrices, and the same measured output y k . It means that ζ k was already computed with (14), thus there is no need to run the Kalman filter twice.

In the algebraic equation ( 22), ζ k is computed through (21) with the Kalman filter ( 14), 19), hence the parameter increment θ is the only unknown, apart from the white Gaussian noise e k . In this linear Gaussian framework, it is well known that the generalized likelihood ratio (GLR) test [START_REF] Basseville | Detection of Abrupt Changes -Theory and Application[END_REF]) for θ = 0 against θ = 0 amounts to

C 0 = C(θ 0 ) is known from the nominal model, Γ k is computed through (23), Φ k is expressed in (
Ω = N k=1 (C 0 Γ k + Φ k ) T Σ -1 (C 0 Γ k + Φ k ) (25a) β = N k=1 (C 0 Γ k + Φ k ) T Σ -1 ζ k (25b) s = β T Ω -1 β. (25c)
The resulting statistics s follows a χ 2 distribution of dim( θ) = n + nr degrees of freedom, central if θ = 0, otherwise non-central with its non-centrality parameter equal to θT Ω θ. The decision for change detection is thus made by comparing s to a positive threshold.

ALGORITHMIC SUMMARY

The computational steps for the derived fault detection method are as follows.

In the nominal reference state:

(1) Obtain the nominal model parameter θ 0 in (7), the associated system matrices A 0 and C 0 in canonical form as in (5a), and the noise covariances Q, R and S in (8). Usually, these quantities are estimated from system identification, e.g., using stochastic subspace identification with the UPC algorithm [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation[END_REF]].

(2) Compute the nominal Kalman gain K and the innovation covariance Σ in (24).

To test if data {y k } k=1,...,N of the current system correspond to the nominal model or not, the test statistics s is computed as follows:

(1) Apply the nominal Kalman filter to the data, obtaining the one-step-ahead state predictions x 0 k in ( 14) and the innovations ζ k in ( 21).

(2) Compute the matrices Ψ k in ( 17) and Φ k in (19), and subsequently Γ k in ( 23).

(3) Compute the test statistics s in (25).

To decide between θ = 0 (no eigenstructure parameter change) and θ = 0, the test statistics s is compared to a threshold. This threshold can be obtained empirically based on realizations of s on data from the nominal state for a given type I error.

APPLICATIONS

Numerical example

Vibration-based damage detection on a simulated massspring chain with eight elements (Fig. 1) is considered as a first application of the derived fault detection method. The matrices M, C and K in (1) of the nominal structural model are defined based on the masses

m 1 = m 3 = m 5 = m 7 = 1, m 2 = m 4 = m 6 = m 8 = 2, stiffnesses k 1 = k 3 = k 5 = k 7 = 1000, k 2 = k 4 = k 6 = k 8 =
500 and a damping ratio of 2% for all modes. Datasets containing output-only time series of accelerations with time step τ = 0.05 s are simulated for different structural states at the four sensor coordinates from white noise excitation at all structural elements. White measurement noise is added with a magnitude of 5% of each generated output signal. 2, where a threshold (red line) is drawn from the nominal state for a 1% type I error. At this type I error, the power of the test for the 5% damage is 97%, and for the 10% damage the power of the test is 100%.

Lab experiment on a beam

Experiments on a PVC beam were carried out by Brüel & Kjaer. The structure's dimensions are 50 cm × 8 cm × 1 cm, being fixed on one side (Fig. 3). Damage is introduced in the beam by drilling small holes. Two damaged states with three holes and with five holes, respectively, are considered.

For both the nominal and the damaged states, acceleration datasets of length N = 295,936 with a sampling frequency of 8192 Hz were recorded under white noise excitation by a shaker. The output data from the nine horizontal sensors on the top of the beam are used in this application. The data are downsampled and decimated by factor 6 to focus on the frequency range of interest, where the biggest damage led to less than 5% decrease in the structure's natural frequencies [START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF]. The increase of the test's non-centrality parameter with a bigger damage is also visible.

CONCLUSIONS

In this work, we have derived a test for the detection of small changes in the eigenstructure of a linear state-space system with application to vibration monitoring of mechanical or civil structures. The test is set up on variables that are easily estimated from system identification and the Kalman filter in the nominal (reference) state of the monitored structure. After transformations of the original detection problem, the standard GLR test is applied in the framework of a simple linear regression subject to additive white Gaussian noises. The normalization of the GLR test taking into account random uncertainties is solely based on the estimated state and output noise covariance matrices, which are standard results of the system identification procedure establishing the nominal model of the monitored structure.

The simple linear regression form of the transformed eigenstructure change detection problem makes possible the implementation of efficient algorithms for on-line monitoring. Moreover, when completed with a finite element model of the monitored structure expressed in a physical parameterization, as e.g. in [START_REF] Döhler | Efficient computation of minmax tests for fault isolation and their application to structural damage localization[END_REF], it is possible to determine the physical parameter(s) responsible for a detected damage. These topics will be the subject of future works.

Fig. 1 .

 1 Fig. 1. Mass-spring chain with four sensors. Stochastic subspace identification with the UPC algorithm [Van Overschee and De Moor, 1996] is used on a simulated

Fig. 3 .

 3 Fig. 3. Experimental setup of the beam.

Fig. 4 .

 4 Fig. 4. Histogram of test values for the beam in the nominal state, and with three and five holes.

  For each structural state, datasets of length N = 10,000 are simulated and the test statistics s is computed. The test values of the nominal state are used to set up an empirical threshold to decide between θ = 0 against θ = 0. The resulting test values for 1000 datasets are shown in the histogram in Fig.

	200
								nominal state
								5% damage
	150							10% damage
	100
	50
	0	0	100	200	300	400	500	600	700
					test value s
	Fig. 2. Histogram of test values for the mass-spring chain
			in the nominal state, and 5% and 10% damage in
			spring 2.
	dataset of length = 200,000 to obtain the nominal
	model parameters. Three different structural states are
	considered for monitoring: the nominal state and two
	faulty (damaged) states with 5% and 10% stiffness de-
	crease in spring 2, respectively.

ACKNOWLEDGMENTS

We thank Dmitri Tcherniak from Brüel & Kjaer for providing the data from the beam experiment.


Maurice Clerc

Stagnation Analysis in Particle Swarm Optimisation or What Happens When Nothing Happens

Introduction

The "classical" PSO version is very simple but the user have to define some parameters (swarm size, neighbourhoods, some coefficients). An adaptive version like TRIBES [START_REF] Clerc | L'optimisation par essaims particulaires[END_REF][START_REF] Onwubolu | TRIBES application to the flow shop scheduling problem[END_REF] does not have this drawback. It is also often more effective when the criterion takes into account only the relation "best value vs number of fitness evaluations". However it is more time consuming, for it performs some intermediate computations in order to take advantage of the information collected during the process. This can be quite annoying, or even unacceptabe for some (quasi) real time applications.

That is why it is still interesting to just improve the classical version, without complicating it too much. The idea here is to mathematically analyse the behaviour of the particles when there is no improvement over several time steps. If the system is seen as a black box that displays only the best positions found so far by the particles, it happens nothing seen from outside, as there certainly happens something inside, which is worthly to study. Better parameters can be derived from such an analysis. Five of the suggested PSO variants are tested on five classical functions, just to bring out some that seems interesting.

Reminder of classical PSO

Principles

Particles are agents moving in a D-dimension search space. The set of this particles is the swarm, whose size is S. At each time step, each particle is described by four features:

-its position -some information links with some other particles (its "informant group", sometimes called its neighbourhood). These "who informs who" links, are often defined once at the beginning, but not topology has been proved to be globally the best one. So we use here a more robust method (see details in [START_REF] Clerc | L'optimisation par essaims particulaires[END_REF]). For each particle K links are randomly drawn at initialisation and then before the next iteration if the previous one didn't improve the best position found by the swarm.

Movement equations for a particle and for each dimension d are given by

( )( ) ( )( )    + ← - + - + ← d d d d d d d d d v x x x g alea c x p alea c v c v 1 , 0 1 , 0 max max 1 (Equ. 1)
where ( )

D g g g ,..., 1 =
is the best position of the best informant.

PSO-0

We need a "reference" algorithm, called here PSO-0. It is the one described above, with the following parameters:

( ) ( ) ( )        ≅ + = ≅ = = = 48 , 1 2 1 72 , 0 2 ln 2 / 1 3 30 2 1 max 1 c c c K S
The first three values has been chosen so that they can be kept in all variants described later. The coefficient max c is derived from 1 c by using the constriction method explained in [START_REF] Clerc | The Particle Swarm-Explosion, Stability, and Convergence in a Multidimensional Complex Space[END_REF] (see Annexe for a quick explanation).

Stagnation phenomena

In classical PSO as defined by equation 1 each dimension is independent. So the analysis can be performed on just one dimension. In order to simplify notations, we write now c instead of max c , and we omit the d index. On the other hand, we will need to explicitely indicate the time step. At last, a given realisation of the random

uniform variable on [ ] c , 0 is now written k c ~,
where k is used to distinguish between realisations. Equation 1 is the rewritten as

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )      + + = + - + - + = + 1 1 1 2 1 1 t v t x t x t x g c t x p c t v c t v (Equ. 2)
We say there is stagnation if p and g (which are now just numerical values) don't change over some time steps. In order to qualitatively observe the behaviour of the particles (along just one dimension, as already noted), we choose some usual values for parameters and initialisations:

If the particle is the (local) best

In order to study the behaviour ot the best particle, let's temporarily set g p = . On figure 1 we can see a typical evolution of the velocity: it tends towards zero, and the particle tends to its best known position, i.e. g. If the particle were intelligent, one could say it thinks "Well, I'm already the best, and the more I see the other don't make any improvement, the less I want to spend my energy to look for a better position". In other words, the particle quite rapidly gives up its search. In practice, it means there are some unuseful fitness evaluations, for anyway the particle (almost) doesn't move anymore.

However the behaviour is depending on 1 c and c. For example, with the values 0.9 and 2.1, we obtain the figure 2, which seems to be a divergent evolution. Of course such a curve can't be a proof, for the velocity may tends towards zero later, but, at least, it is a strong presumption. The lesson that can be drawn is there probably exists an intermediate parameter set so that the particle keeps some activity but not diverges.

Behaviour of the (local) best particle in case of stagnation, with some unusual parameter values

If the particle is not the best one

Actually, the examination of what happens with a non best particle, i.e. with g p ≠ , will suggests some possible solutions. The behaviour is indeed very different, even with usual parameter values. As we can see on figure 3 the velocity is constantly modified , without tends to zero and the position is permanently oscillating around a "mean" one.

A deeper examination shows that the module of the velocity follows a statistical law that is exponentially decreasing (cf. figure 4): small oscillations are frequent, more important ones are more rare, big ones are even more rare, etc. In other words, the longer the stagnation, the more probably big moves will occur. This is a quite desirable behaviour for then, even during stagnation, it performs a kind of local search that is usually fruitful: the particle will probably find a better position. So a quite obvious (small) improvement could be to do something so that the best particle has the same behaviour.

Figure 3. PSO-0. Behaviour of a non best particle in case of stagnation, with some usual parameter values. In case of long stagnation (over 100000 time steps here) a non best particles oscillates according to an exponential decreasing probability distribution: small moves are quite frequent, big ones are rare but do exist

Stagnation analysis

.

General iterative equation and interpretation

In the particular case 0 = c , we obviously have ( ) ( )

t v c t v 1 1 = +
, and the convergence criterion is simply 1 1 < c (we assume it is non negative). In the general case, from equation 2 we can successively write

( ) ( ) ( ) ( ) ( ) (
)

1 1 2 4 3 1 + - + + - + = + t x g c t x p c t v c t v ( ) ( ) ( ) t x t v t x + + - = + - 1 1 ( ) ( ) ( ) ( ) g c p c t v c t v c c t x 2 1 1 2 1 1 1 - - - + + = - ( ) ( ) ( ) ( ) ( )w g p t v q c t v q s c t v 1 2 1 1 - - - + + - = + (Equ. 4) with        - = + + = + = q c c w c c c c q c c s 1 3 2 1 4 3 4 3 (Equ. 5) Quantities s ~, q ~, et w ~ are realisations of the random variables S, Q et W. Let's also define Q S c Z + - = 1
. The general iterative equation is then based on the three random variables Z, Q et W . It can be rewritten in a summarized form, and by shifting the iteration ranks:

( ) ( ) ( ) ( )W g p t Qv c t Zv t v - + - - = + 1 1 1 (Equ. 6)
Probability densities of S and Q can quite easily be formulated (cf. Annexe). The ones of Z and W are extremely complicated, but anyway we can have an idea of their "shape" (by simulations), and their main features (minimum, maximum, mean) can be computed (cf. Table 1. Details and figures are in Annexe). Note that Z and Q have an infinite support, with a density tending towards zero for high values, and that explains the exponential decline that we have seen in figure 3.

. Minimun Maximum Mean Z c c 2 1 1 - + ∞ ( ) 2 ln 2 1 + -c c Q 0 ∞ ( ) 2 ln 2 W c - c 0 (std. dev. ( ) ( ) 12
2 ln 14 11c )

Table 1. Main features of the three random variables of the general iterative equation

Let us try now to give the "meaning" of the three terms. The last one (W) can be seen as a "noise". The bigger g pthe more important the noise. This is intuitively quite sastifying: the more the particle is far from the best position it knows, the more it tends to try big moves. Note that there is a finite maximum possible value (and a finite minimum one). It might be better to replace this noise by a classical Gaussian one.

As

Q c 1 -
is always negative and as Q is parameter free, this term can be seen as a "back force", whose intensity is only depending on 1 c .

At last Z, the first one, is usually sometimes negative (back force) sometimes positive (forth force). However we can note that the back force may exist only if c c 2 1 1 -+ is negative. Again this suggests a possible relationship, that is to say

2 1 1 + = c c (Equ. 7)
This would clarify the role of each random variable: noise for W, back force for Q, and only forth force for Z. However, in order to find some other constraints on parameters we have to go further in our analysis.

If the particle is the best one

For the best particle we have g p = and we clearly see from equation 6 that there is then no noise component. In passing, it explains why the particle may converge for some parameter values. If we consider two consecutive time steps (don't forget we are in stagnation, so p is constant), the equation can be rewritten in a matrix form

( ) ( ) ( ) ( ) ( ) ( )       - =       -         - =       + 1 1 0 1 1 1 t v t v M t v t v q c z t v t v t t t (Equ. 8)
Then the behaviour is defined by

∏ = t k k M 0 ~.
To study it is quite difficult (see Annexe), but the final result is a relation between [START_REF] Clerc | L'optimisation par essaims particulaires[END_REF] c and c that we can see on figure 4. Any pair of parameters in the convergence domain is usable. However, in iterative optimisation, it is often a good thing to work "on the edge of chaos". So points on the frontier curve should be better. Moreover, in order to enlarge the local search, 1 c and c should be "both as big as possible". There are several way to formulate this, for example by maximizing the quantity c c 1 . Finally, the most promising parameters are on the curve and "around" the point

( ) ( )      ≅ + = ≅ = 39 . 2 2 ln 2 1 72 . 0 2 ln 2 1 1 c c ( Equ. 9)
Remember it is only for stagnation phases. If we want to use just one set of parameters, no matter there is stagnation or not, it may be a good idea to choose say 1 c according to the equation 9, and c according to equation 7, for example (see variant 3PD-PSO-0 below).

Figure 4. In case of stagnation best parameters are on the frontier of the convergence domain and, on the other hand, both "as big as possible"

If the particle is not the best one

Let's remind that our hypothesis is stagnation: neither p nor g are modified. Then, when the particle is not the best one, that is to say when g p ≠ , it can't converge, because of the "noise" component W. It is indeed what we observe, no matter how small 1 c is. Actually, even for 0 1 = c the particle is oscillating. However it can diverge it this coefficient is too big.

More precisely, for the mean of W is zero, the conditions that we have found for the best particle are still valid, except that they don't imply convergence but non divergent oscillations. This is quite a good thing for this behaviour performs a kind of local search..

Generalisation and variants

The three distributions

We can now rewrite equation 6 in a more general way

( ) ( ) ( ) ( ) 3 2 1 1 1 R g p t v R t v R t v - + - - = + (Equ. 10)
where the three probability distributions have the following features 1 R is positive unimodal on a finite support 2 R is positive unimodal on a finite support 3 R is symmetrical and unimodal with a mean equal to zero

We have seen that we can choose

W R = 3 , Q c R 1 2 =
, and Z R = 1 , with the condition

( ) 2 / 1 1 + ≤ c c
, but some others are clearly possible. In particular, as 1 R is for the "forth force" and 2 R for the "back force", we can manipulate c and 1 c so that they are quite similar. We will see a simple example below. Starting from PSO-0, we will now progressively add just some small modifications. Don't forget that although we are looking for improvements we also want to keep the algorithm as fast as possible.

There will be two kinds of modification that takes advantage of our stagnation analysis:

-the ones that just define some new parameter sets, that are used no matter there is no stagnation or not (see below variants 0, 0' and 1)

-the ones that add a rule like "if there is stagnation, then ..." (see below variants 2 and 3)

Theoretically, stagnation occurs as soon as p and gdon't change for two time steps. However, as we have seen, it is not a good idea to do something special immediately, at least as long as the best particle is moving "enough". So we have first to define how long should be the stagnation before to modify the moving strategy.

Defining a stagnation threshold

The underlying idea is to say there is indeed time do try something different if:

the local best particle has been enough time to be informed (either directly or undirectly) by a possible better particle

-it has nevertheless not improved its best known position Note that this criterion is not pertinent for the global best, but in order to keep the algorithm simple, we neglect this particular case.

The probability that a given particle A is still not informed after t time steps that there exists a better particle G is given by ( )

Kt t S A G p       - = → ¬ 1 1 (Equ. 11)
where S is the swarm size and K the number of information links drawn at random (at each time step, for there is no improvement. For detailed formulas of "temporal connexity", see [START_REF] Clerc | L'optimisation par essaims particulaires[END_REF]) . If we want that this probability is smaller than a threshold ε the number of time steps should be at least equal to , we would now say there is stagnation if the local best particle didn't improve its position from at least 91 time steps. This is the value used in some variants below. The additional parameter ε is unfortunately empirical: this is the price to pay to be sure the algorithme remains fast.

Test bed

We have now in hand all what we need to test some variants on a some functions, in order to compare them. The five functions used here are very classical, and the minimum is always zero. However, the dimension (30), the search spaces and the required accuracies are chosen so that our "reference" algorithm PSO-0 has a success rate never equal to zero nor to 100% (over 100 runs with at most 40000 fitness evaluations). This success rate can therefore be chosen as a comparison criterion, for each function. For a quick comparison the mean success rate is even enough. In case of a tile (or almost), we will prefer the variant that gives the smaller standard deviation, for it is then probably more robust.

Function

Search

Six variants

3PD-PSO is for "Three Probability Distributions Particle Swarm Optimisation". The following variants are just a sampling of the most obvious ones suggested by stagnation analysis. In particular we don't modify at all distribution 2 R , which is still equal to Q c 1 , and 1 R is also equal to Z, just possibly translated. Table 3 summarizes the results.

Function 3PD-PSO-0 3PD-PSO-0' 3PD-PSO-1 3PD-PSO-2 3PD-PSO-

3PD-PSO-0

This variant is simply PSO-0 with some different coefficients. That is to say stagnation is not checked. So we need the find a "compromise", some parameters that are not bad if there is stagnation, and not bad if there is no stagnation. As already suggested, according to the above analysis, a promising choice is

( )        ≅ + = ≅ = 86 . 0 2 1 72 . 0 2 ln 2 1 1 1 c c c
so that the random variable Z (the "forth force") never gives negative values. The movement equation is then directly derived from equation 1 and can be written for dimension d

( )( ) ( )( ) d d d d d d x g alea x p alea v v - + - + ← 1 , 0 86 . 0 1 , 0 86 . 0 72 . 0 (Equ. 13)
Other parameters are the same (30 particles, 3 informations drawn at random for each particles). Effectiveness is significantly improved.

3PD-PSO-0'

Here we keep

( ) 72 , 0 2 ln 2 1 1 ≅ = c
, but we choose c so that the mean of Z (for the "forth force") is equal to the one of Q c 1 (for the "back force"). It gives

( ) ( ) 108 . 1 1 2 ln 2 2 ln 2 1 ≅ - + = c
. The velocity is then updated by

( )( ) ( )( ) d d d d d d x g alea x p alea v v - + - + ← 1 , 0 108 . 01 1 , 0 108 . 1 72 . 0 (Equ. 14)
Effectiveness is more improved than with 3PD-PSO-0, and the standard deviation is also reduced.

3PD-PSO-1

For

( ) 72 , 0 2 ln 2 1 1 ≅ = c
, stagnation analysis suggested a possible c (i.e. max c in equation 1) equal to 0.86. The idea is here to define another plausible value, and to choose c at random between them for each move. For example value 1.48 used in PSO-0 is not arbitrary: it comes from a deterministic convergence analysis (cf. Annexe). We can then try the variation interval U being two uniform random variables on [ ] 1 , 0 . This distribution is obviously not uniform anymore, as we can see on figure 10. Intuitively it is more satisfying, for big moves are now less probable. The improvement is equivalent to the one given by 3PD-PSO-0 , although the standard deviation is slightly higher.

[ ] ( )         + + = 2 1 , 2 1 ,

3PD-PSO-2

In this variant, we define two sets of parameters, one for "normal" case, and one for stagnation. In the first case, we use the parameter set of 3PD-PSO-0', and in the second case, we choose the "central point" of the promising area that we have seen on figure 4:

( ) ( ) ( )        ≅ - + = ≅ = 108 .

3PD-PSO-2'

We have seen that in case of stagnation the local best particle has a special behaviour that rapidly induces some unuseful fitness evaluations. Using two set of coefficients is already a way to change this, but can also modify the movement equation so that the "noise" term is not null.

Let G be the local best particle. In order to compute its new velocity, we can consider all best known positions in the swarm (in short, all "p"), and to keep the ones that are better that the one of G (i.e. g). Then we choose one at random, say 2 g , and in the movement equation we replace g by 2 g . So, except for the global best, the movement equation is now similar to the one of a non best particle.

By adding this modification to 3PD-PSO-2, the effectiveness it more or less the same and the standard deviation is a bit smaller.

3PD-PSO-3

Finally, we have seen that the "noise" W is vaguely similar to a truncated Gaussian curve. Therefore, in case of stagnation, it is tempting to replace it by a true Gaussian distribution. In order to not have to define a new parameter, we can choose this Gaussian distribution with the same standard deviation that W, i.e.

( ) (

)

12 2 ln 14 11 - = c σ (almost 3 c ).
The algorithm is then a bit more complicated. We need to memorize velocities at time t and t-1, and also the random values that are used. Then, in case of stagnation velocity is directly computed by equation 6. In practice, it gives

( ) ( ) ( ) ( ) ( ) ( ) σ , 0 _ 1 1 1 1 normal alea g p t v q c t v q s c t v d d - - - - + - = + (Equ. 18)
The process is not really slowed down. However applying this method to 3PD-PSO-0' doesn't improve it (the result is the same if applied to 3PD-PSO-2).

Conclusion

There is a slight advantage to variant 3PD-PSO-2' (best mean and seems a bit more robust). Let's summarize the differences compared to PSO-0:

-two sets of coefficients, one when everything goes well, and one where there is stagnation -for the local best the set of informants is the whole swarm, and if possible a better than itself it chosen (at random in case of a tile)

Other variants based on stagnation analysis are of course possible. In particular we could more completely manipulate the "forth force" and "back force" components, for example by using Beta distributions. And more tests are needed to choose the best variant.

Annexe

Mathematical developments below have been simplified. In particular the details of some tedious probability densities are not given, for they are not really more helpul than their graphical representations to understand this study.

Relationship between coefficients according to deterministic analysis

In [START_REF] Clerc | The Particle Swarm-Explosion, Stability, and Convergence in a Multidimensional Complex Space[END_REF] several constriction methods are defined, that guarantee the non-divergence of the swarm. For historical reasons, the most commonly used gives the two coefficient 1 c and max c as depending on an unique parameter ϕ :

       = - + - = 2

Figure 1 .Figure 2 .

 12 Figure 1. PSO-0. Behaviour of the (local) best particle in case of stagnation, with some usual parameter values

  In practice, to compute the new velocity of a particle, max c is first chosen at random each dimension d equation 1 is applied. Formally, it is equivalent to apply equation 21 realisations of the random variable B, defined by

  This variant is again slightly better both for the mean and the standard deviation.

and 2 X

 2 be two uniform random variables on [ ] 0 else. One can distinguish two cases (we do suppose here c is not null): mean is the sum of the means of 1 X and 2 X , i.e. c S = ˆ. Actually S can be rewrittencU S =where U is a uniform random variable on [ ]

2 S=

 2 Figure 5. S density for c=1.48

Table 2 . Test bed and success rate with PSO-0 (100 runs for each function, and at most 40000 fitness evaluations for each run)

 2

	Parabola (Sphere) Griewank Rosenbrock Rastrigin	space 20 , 20 ] 30 300 , 300 -[ -[ ] 30 [ 10 , 10 -] 30 [ ] 30 12 , 5 , 12 , 5 -	Accuracy 10 -0.0001 25 35	Success rate with PSO-0 32% 44% 21% 12%

Table 3 . Results with some variants. When the success rate is equal to 100%, the mean number of fitness evaluations is about 18000

 3

	2' 3PD-PSO-3

Acknowledgements

Jean-Pierre Levrel (newsgroup fr.sci.maths) Ray Koopman (newsgroup alt.sci.math.probability)

Note that Q is parameter free. The mean can be computed by noting that ( )

, as a product of two independent random variables. This mean is then the one of ( )

. The density of this last random variable is

Its mean is then

W study

The random variable W can be rewritten

So we immediately see that it is symmetrical and that its mean is zero. In order to compute the maximum value, one can note that it is reached when 1 X is null and that in this case W is equal to 3 X . This maximum is then equal to c. The formulas that give the density are quite complicated, for they make use of the order 2 Jonquières

. It is enough here the visualize it and to compute the standard deviation

Figure 7. W density for c=1.48, and best Gaussian approximation

Z study

The random variable

. We have seen that the mean of ( )

Therefore the mean of Z is ( )

, but this formula couldn't have been directly written, for random variables S and Q are not independent. Here again, we just visualize the density curve, without giving the too complicated formulas.

B study

The random variable B can be rewritten

The probability density is given by

Convergence of an iterative random sequence

What follows is just a sketch of the proof.

A and B are two random variables. ( ) ( ) ( )

This can be rewritten

We want to study the convergence of y to zero. Let's define

, a realisation of the random matrix

. So the condition is simply 0 ˆ< G . We have

In the particular case of stagnation analysis, we have

Note that there is trick here, for k m ~and 1 ~+ k m are not independent. However k m ~and 2 ~+ k m do are, so the condition is still valid. Unfortunately the result is quite complicated. A good quadratic approximation can be found by saying that the frontier of this convergence domain contains the three points ( ) ( )

, and ( ) 0 , 1 (see figure 4), and by empirically estimating two parameters α and β :

( )

. So, as soon as we have

the system begins to converge. That is why the value ( )

c is smaller than 1 χ , convergence begins immediately. If it it bigger than 1 χ (but still smaller than 1), convergence may begin only after ( ) ( )

Stagnation in FIPS (Fully Informed Particle Swarm)

When more than just two particles inform a given one, the motion equations can be written

where 1 ~k c is a realisation of a random uniform variable on [ ] c , 0 .

If no k

p is modified , we can deduce the formulas are on the whole quite similar. The formula for ω ~ is of course more complicated, but we can note that if all k p are equal then we still have 0 ~= ω . In particular it would be true for the global best: the "noise" component would again disappear. However, in FIPS each particle is not informed by itself, and we see here why the probability there is no noise component is then almost null.


Vincent Chabannes

Alexandre Ancel

Julien Jomier

Christophe Prud'homme

An Open Platform to reconstruct vessels from MRI images and simulate blood flows to ultimately provide Virtual Angiographies

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

AN GI OTK

Fr om m edical im ages to num er ical sim ulations

ANGIOTK

The platform is still a prototype, but reached a level of maturity that requires very little human intervention. The overall objective is to make available, to the medical community, tools for computer aided modeling for conducting experiments in silico, otherwise difficult or impossible to carry on patients. Feel++ is an open-source C++ library dedicated to solve partial differential equations using a language very close to the mathematics embedded in C++.

Its features include: -arbitrary order Galerkin methods from 1D to 3D -various finite elements H, L, H(div), H(curl) -seamless parallel computing -scale up to tens of thousands of cores -In-Situ Visualisation

A

  specific aim is to create an ecosystem around the platform composed of academic, research centers and enterprises and to provide services such as training, specific developments, consulting or deployment The platform is operated by Cemosis and Kitware. An Open Pla tform to reconstruct vessels from MRI ima ges a nd simula te blood flows to ultima tely provide Virtua l Angiogra phies V. Chabannes, A. Ancel, J. Jomier, C. Prud'Homme Rencontres INRIA-Industries Santé, Bordeaux , Oct. 13 2015Cont act s:christophe.prudhomme@cemosis.fr julien.jomier@kitware.com

  We have developed the open-source platform AngioTK that implements the complete pipeline using the following open-source software Imaging: VTK & ITK http:/ / www.itk.org Filtering: RORPO http:/ / path-openings.github.io/ RORPO/ Extraction: VMTK http:/ / www.vmtk.org/ Mesh Generation: GMSH http:/ / geuz.org/ gmsh Blood Simulations: Feel ++ http:/ / www.feelpp.org MRI Simulations: JEMRIS http:/ / www.jemris.org A Softwa re Pipeline Devel opment https:/ / github.com/ vivabrain/ angiotk Forum https:/ / gitter.im/ vivabrain/ angiotk Open-Source License BSD

Funding Consortium 1 2 Websit e http:/ / www.feelpp.org Devel opment https:/ / github.com/ feelpp/ feelpp Forum https:/ / gitter.im/ feelpp/ feelpp


Luca Maggiani

Cedric Bourrasset

Matteo Petracca

Francois Berry

Paolo Pagano

Claudio Salvadori

HOG-Dot: a Parallel Kernel-Based Gradient Extraction for Embedded Image Processing

In this paper we propose HOG-Dot, a method for the direct computation of the polar image gradients coordinates from the pixels values. The proposed algorithm, to be used as the first step of the Histogram of Oriented Gradient (HOG) pipeline, approximates the exact gradient with its projection onto a versor chosen among the projection plane set. Instead of non linear computations, the HOG-Dot method exploits linear operations while introducing a bounded approximation error with respect to other HOG approaches, thus resulting a more suitable solution for embedded devices. Concerning the state of the art, it also achieves improved accuracy with the mathematical spatial gradient formulation.

I. INTRODUCTION

Embedded image processing devices, namely Smart Cameras (SCs), are nowadays considered as a powerful solution to extract complex information from the surrounding environment in which they are deployed. In designing SC one of the biggest efforts is the porting of complex and computational intensive computer vision pipelines to resource constrained embedded devices able to reduce deployment costs. In such a context, the Histogram of Oriented Gradient (HOG) algorithm [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] represents the state of the art for the pedestrian detections, allowing the classification of images windows represented by using histograms of local spatial derivatives through machine learning algorithms (i.e., Support Vector Machine [START_REF] Cortes | Support-vector networks[END_REF]). The HOG pipeline is composed by a gradient extraction followed by a local spatial aggregation histogram. The resulting histogram is then normalised and compared to a model reference within the classifier. In particular, the spatial gradient extraction is based on intensive pixel-wise operations usually exploited with non-linear operators [START_REF] Bhardwaj | A survey on various edge detector techniques[END_REF]. Thus its deployment into resource constrained embedded devices does not guarantee real-time processing performance. Within the HOG algorithm, spatial gradient extraction is usually computed with square root and arctangent operators from the 1-D spatial derivatives. In Eq. 1 the 1-D components are computed as:

G x = ∂I ∂x = I x + 1, y -I x -1, y (1)
G y = ∂I ∂y = I x, y + 1 -I x, y -1
The gradient magnitude ∇I(x, y) and orientation angle θ are then computed for each pixel as follows:

∇I(x, y) = G 2 x + G 2 y (2) θ = arctan G y G x
The resulting gradient orientation θ is then sampled over a set of discrete samples, called bin, to create a dense grid of uniformly spaced histogram cells.

Since gradient extraction represents a common image processing technique, hardware and software implementations, especially for embedded devices, have been proposed in literature. In [START_REF] Katic | A relative imaging cmos image sensor for high dynamic range and high frame-rate machine vision imaging applications[END_REF] several VLSI gradient extraction implementations have been reviewed. Such solutions deploy gradient magnitudes evaluations directly in a custom CMOS camera focal plane as pre-processing step. In [START_REF] Kadota | Hardware architecture for hog feature extraction[END_REF] a Look-Up table (LUT) implementation has been proposed. It exploits a set of precomputed Look-Up Tables (LUTs) to approximate both the square-root and the arc-tangent. The results are used to index a 2D LUT and to address the preloaded magnitude and angle values. In such implementation the memory footprint grows accordingly to the size of the 2D LUT, and consequently to the gradient resolution. More recently, in [START_REF] Lee | Accelerating histograms of oriented gradients descriptor extraction for pedestrian recognition[END_REF] the polar conversion is computed implementing both the arc-tangent and the square-root functions with iterative hardware operations, thus reducing significantly the data throughput.

In this paper we present HOG-Dot, a kernel-based gradient extraction algorithm. It permits to compute directly the gradient polar components by deploying linear operations. Instead of evaluating an approximation of the Eq. 2 or deploying a custom circuitry as shown in [START_REF] Katic | A relative imaging cmos image sensor for high dynamic range and high frame-rate machine vision imaging applications[END_REF], we rather propose a mathematical reformulation which offers a bounded approximation error with respect to the exact gradient. By modifying the number of bins, the resulting error is then controllable and deterministic.

The rest of this letter is organized as follows. First, in Sec. II the proposed spatial gradient expression is mathematically evaluated. In Sec. III our method is then compared with other state of the art gradient extraction formulations and finally Sec. IV concludes the paper.

II. METHODOLOGY

One of the main issue in developing the HOG feature extraction for embedded SC is the polar coordinates computations through non linear floating-point operations. Even though floating-point arithmetic could be exploited even in resource constrained devices, as shown in [START_REF] Salvadori | An optimisation of gaussian mixture models for integer processing units[END_REF] its implementation requires a non-negligible amount of time thus decreasing the output throughput. Thus, in this paper we propose a linear method capable of approximating the square-root and the arctangent. It is based on a set of linear operations that can be easily implemented within embedded device. Let us consider a generic vector v, its module can be seen as a dot product with its parallel versor:

v = v • îθ ( 3
)
where θ is the angle and then the direction of the versor parallel to v. By exploiting this property, we can state that the exact gradient is equal to:

∇I(x, y) = ∇I(x, y) • îθ ( 4
)
where θ is the direction of the vector ∇I(x, y). Moreover, considering a generic versor îk , the following relation is valid:

∂I ∂ îk (x, y) = ∇I(x, y) • îk (5)
Since the HOG algorithm considers a discrete angle domain that consists of equidistant samples within the interval [0, π), we consider a limited set of N directions, represented by the following of versors (illustrated in Fig. 1):

B î = { îk (θ k ) | θ k = kπ/N, k = 0, . . . , N -1} (6)
where each versor can be represented in cartesian coordinates as:

îk = x cos θ k + ŷ sin θ k (7) î0 îk îk+1 îN-1 θ k
Fig. 1. The versors beam used to compute the dot product.

According to the above reported equations, the intensity of the gradient projection ∇I(x, y) onto each versor of the set B î can be evaluated as follows:

∇I(x, y) • îk = ∂I ∂ îk (x, y) = cos θ k ∂I ∂x + sin θ k ∂I ∂y (8)
Since we consider the HOG case, the partial derivative functions are expressed in terms of pixel intensity and position as in Eq. 1.

Thus, by merging Eq. 8 and 1, the value of ∂I ∂ îk is then derived as follows:

∂I ∂ îk (x, y) = cos θ k I(x + 1, y) -I(x -1, y) + sin θ k I(x, y + 1) -I(x, y -1) (9)
As in Eq. 9, the ∂I ∂ îk value represents the magnitude of the gradient projection onto the versor îk . It can be computed by convolving an image with a kernel matrix filled with constant elements (θ k are a-priori defined as well as the cos θ k and sin θ k values).

Among the k computed projections, the argmax function is used to extract the greatest gradient projection as results of magnitude comparisons, as shown in Eq. 10. The extracted k projection also represents the closest gradient approximation provided by our method.

k = argmax ∂I ∂ îk (x, y)|k = 0, ..., N -1 (10)
as results of Eq. 10, the angle θ of the ∇I(x, y) vector belongs to the following interval:

θ ∈ θk - π 2N , θk + π 2N (11)
Thus, the HOG-Dot method approximates the exact gradient ∇I(x, y) as follows:

∇I Dot (x, y) = ∂I ∂ îk (x, y) , θk (12)
where θk = kπ/N . According to the original HOG implementation [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] ∇I HOG (x, y) is expressed as:

∇I HOG (x, y) = ( ∇I(x, y) , θk) (13)
Considering now Eq. 12 and Eq. 13, the vector obtained by our method is overlaying the HOG one but their magnitude differs. The Dot gradient is indeed representing the projection of the exact gradient over the îk versor (as in Eq. 14) while the HOG version has the same magnitude as the exact gradient. In Fig. 2 the geometric correspondences between both methods are shown. As result of Eq. 4 and Eq. 12, the formal relation between ∇I HOG (x, y) and ∇I Dot (x, y) modules is then expressed in Eq. 14. In particular, the ∇I Dot (x, y) represents a round down approximation of the ∇I HOG (x, y) gradient.

θk + π 2N θk -π 2N îk îk +1 îk -1 ∇ I ( x , y ) ∇ I D o t ( x , y ) ∇ I H O G ( x ,
∇I Dot (x, y) = ∇I(x, y) cos(θ -θk) (14)
In the following Section, the proposed approximation is quantitatively evaluated with respect to other HOG implementations.

III. ERROR EVALUATION

In this section the HOG-Dot gradient is first compared with the original HOG algorithm (Sec. III-A). In Sec. III-B an extended comparison is then shown with respect to the exact gradient formulation. As assessment metric we selected the Error Vector Magnitude (EVM) [START_REF] Michael | Evm calculation for broadband modulated signals[END_REF], an index commonly used in telecommunications to measure the difference between transmitted and received constellation points (distance between reference and real points). In mathematical terms EVM can be expressed as:

EV M = P error P ref ( 15
)
where P error is the power of the error vector and P ref is the power of the reference measure. For our purposes the EVM can be rewritten as:

EV M = e 2 ∇I ref (x, y) 2 = e ∇I ref (x, y) (16)
where e is ∇I ref (x, y) -∇I real (x, y) , the difference vector between the reference and its approximated version. In the following of the section, the notation EV M A B is used to specify the reference (A) and the real (B) point considered in the assessment, e.g., EV M HOG Dot is evaluated between the HOG and the HOG-Dot approximations.

A. HOG-Dot versus HOG

As in Eq. 14 the difference vector between ∇I HOG (x, y) and ∇I Dot (x, y) is straightforward evaluated.

e HOG Dot = ∇I HOG (x, y) -∇I Dot (x, y) = = ∇I(x, y) 1 -cos(θ -θk) (17)
Thus, the EV M HOG Dot is a deterministic function of θ -θk, and, consequently of N (θk = kπ/N ). The EV M HOG Dot is then described as follows:

EV M HOG Dot = 1 -cos(θ -( kπ/N )) (18)
The EV M HOG Dot function is then bounded within the interval 0, 1 -cos( π 2N ) . It reaches the lower bound value when θ = θk and the maximum value when θ = π 2N (the border of the interval). Moreover, since EV M HOG Dot is function of N , the maximum approximation error can be controlled by fixing an appropriate angle resolution as in Eq. 11. In Fig. 3 To complete the comparison, maximum and average errors, respectively annotated max(EV M HOG Dot ) and avg(EV M HOG Dot ), are reported as function of N in Fig. 4. By increasing the angle sampling parameter N , the difference between θ and θk decreases, thus the maximum EV M decreases as well. Given all values on its period, the average EV M is evaluated as follows:

avg(EV M HOG Dot ) = 2 T T /2 0 EV M HOG Dot dθ (19) = 1 - 2N π sin( π 2N )
where T is equivalent to the EV M HOG Dot period π N . The average value is clearly bounded by Eq. 18 and decreases when N increases. Although we exploit the gradient extraction with a linear kernel based method, our results are comparable with those obtained by the HOG technique. As shown in Fig. 4, by fixing N equal to 8 as Dalal et al. suggests, the resulting average EVM is close to 0.5%. Whether an improved accuracy is requested, the sampling factor N can be further increased to match the output specifications. The higher is N , the lower EV M HOG Dot results. As far as N increases, the ∇I projections are closer to the reference ∇I HOG vector, which asymptotically will coincide.

B. HOG methods versus ∇I

In Sec. III-A, the HOG-Dot technique has been compared with the original HOG gradient extraction. In this Section HOG-Dot and other state of the art gradient extraction techniques are compared with the exact gradient formulation ∇I through the EVM metric. To better understand the main differences between the the HOG-Dot and HOG methods, ∇I Dot and ∇I HOG are reported in Fig. 5, in which the exact gradient ∇I is also shown. In order to compare the considered methods with respect to the ∇I, in this occasion the EVM has been evaluated with respect to the exact ∇I. Thus, according to Eq. 12, ∇I Dot represents the projection of ∇I along the versor îk and ∇I HOG overlies the îk versor with the same ∇I module. The resulting EVM can be expressed as follows:

θ -θ k ∇ IHO G ∇ I ∇ IDo t îk
EV M ∇I Dot = sin(θ -θk) (20)
where EV M ∇I Dot is the magnitude of the difference vector between ∇I and ∇I Dot . As shown in Fig. 5, the EV M ∇I ∇I HOG is then:

EV M ∇I HOG = 1 -cos(θ -θk) 2 + sin(θ -θk) 2
(21) Thus, by applying Eq. 20 to Eq. 21 the final EV M ∇I HOG is then obtained.

EV M ∇I HOG = 1 -cos(θ -θk) 2 + EV M ∇I Dot 2 (22)
According to Eq. 21, the original HOG gradient approximation proposed by Dalal et al. presents an EV M ∇I HOG higher than EV M ∇I Dot . In Lee et al. [START_REF] Lee | Accelerating histograms of oriented gradients descriptor extraction for pedestrian recognition[END_REF] a gradient approximation has been implemented within an hardware circuitry. In this case the EV M ∇I Lee can be expressed as:

EV M ∇I Lee = 1 ∇I ∇I - G x + G y + √ 2G max 1 + √ 2 (23
) where G x and G y are the gradient 1-D components as in Eq. 1 and G max is the greater absolute value among them. Finally, the Bhardwaj and al. [START_REF] Bhardwaj | A survey on various edge detector techniques[END_REF] implementation is considered. In this case, the gradient magnitude is approximated by directly adding the G x and G y components. The EV M ∇I Bhardwaj is then straightforward evaluated as follows:

EV M ∇I Bhardwaj = cosθ k + sinθ k -1 (24)
where θ k is the HOG binning angle as in Eq. 6.

The proposed algorithm is numerically evaluated with respect to the other methods. Simulations have been performed over different datasets (INRIA [START_REF]Person dataset[END_REF], CVC [START_REF] Cvc | CVC-01 Pedestrian Dataset[END_REF]) by calculating the spatial gradient with the above presented techniques. The maximum, the average and the variance of the EVM, respectively reported as max(EV M ), avg(EV M ) and σ 2 (EV M ) in Table I, have been evaluated for each method. The comparison with the exact gradient shows a minimum EVM error of 19,5% mainly due to the HOG binning operation which introduces an angular error. With this respect, our proposed solution outperforms the other state of the art techniques. As a matter of comparison, the evaluation has been performed with N equal to 8 for all the considered methods. Nonetheless our method involves only linear computation, the resulting EV M ∇I Dot is smaller than others with respect to the ideal gradient. With respect to the method proposed by Lee et al., HOG-Dot achieves better results in terms of the maximum EVM and a reduced average error. With respect to the Bhardwaj et al. method, the HOG-Dot performs significantly better and reduces the maximum EVM by 25%. Finally, the comparison between the Dalal et al. and our technique essentially provides the same EVM performance, with the limited improvement as expected by Eq. 22. Although the proposed linear approximations, HOG-Dot shows better results with respect to the other state of the art techniques. By exploiting only linear computations, HOG-Dot is then suitable for constrained embedded devices, where non-linear, square root and arctan methods will be reducing the processing performance.

IV. CONCLUSIONS

In the paper the HOG-Dot technique has been proposed. The HOG-Dot is a kernel-based gradient extraction technique capable to compute directly the polar coordinates of image gradients by linear approximating the exact gradient. As comparison with the other state of the art algorithms, HOG-Dot shows a bounded and controllable approximation error. Even though our method introduces an error with respect to the exact mathematical formulation, it shows good accuracy results, outperforming the other state of the art HOG methods. Moreover, due to its internal parallel architecture, HOG-Dot is suitable for hardware oriented implementations. As future work the HOG-Dot will be implemented in a FPGA-based SC to evaluate the latency performance and the impact on the approximation error within the HOG pipeline.

Fig. 2 .

 2 Fig. 2. The dot-product between ∇I(x, y) and the "closest" îk , compared with Dalal et al. technique.

Fig. 3 .

 3 Fig. 3. The EV M HOG Dot evaluation as function of θ -θk.

Fig. 4 .

 4 Fig. 4. Maximum and average EV M HOG Dot as function of N .

Fig. 5 .

 5 Fig. 5. ∇I, ∇I HOG and ∇I Dot comparisons.

ACKNOWLEDGMENT

This work has been sponsored by the French government research programme "Investissements d'avenir" through the IMobS3 Laboratory of Excellence (ANR-10-LABX-16-01), by the European Union through the program Regional competitiveness and employment 2007-2013 (ERDF Auvergne region), and by the Auvergne region.