Agustín Jiménez

Basil Mohammed Al-Hadithi

Juan Pérez-Oria

Luciano Alonso

Optimal Control Using Feedback Linearization for a Generalized T-S Model

In this paper, a fuzzy feedback linearization is used to control nonlinear systems described by Takagi-Suengo (T-S) fuzzy systems. In this work, an optimal controller is designed using the linear quadratic regulator (LQR). The well known weighting parameters approach is applied to optimize local and global approximation and modelling capability of T-S fuzzy model to improve the choice of the performance index and minimize it. The approach used here can be considered as a generalized version of T-S method. Simulation results indicate the potential, simplicity and generality of the estimation method and the robustness of the proposed optimal LQR algorithm.

Introduction

Feedback linearization has been used successfully to address some practical problems. These include the control of helicopters, high performance aircraft, industrial robots and biomedical devices. More applications of the methodology are being developed in industry [START_REF] Slotine | Applied nonlinear control[END_REF].

It is well known that a robust a effective controller requires an accurate and efficient model. In [START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF], an interesting method is developed to identify nonlinear systems using input-output data. The main problem encountered is that T-S identification method can not be applied when the MFs are overlapped by pairs. Nonlinear control systems based on the T-S model have attracted lots of attention during the last twenty years (e.g., see [START_REF] Barragán | A general methodology for online TS fuzzy modeling by the extended Kalman filter[END_REF], [START_REF] Gang | A survey on analysis and design of model-based fuzzy control systems[END_REF] and [START_REF] Hseng | Fuzzy Bilinear Model and Fuzzy Controller Design for a Class of Nonlinear Systems[END_REF]). It provides a powerful solution for development of function approximation, systematic techniques to stability and design of fuzzy control systems in view of fruitful conventional control theory. They also allow relatively easy application of powerful learning techniques for their identification from data.

T-S fuzzy models are proved to be universal function approximators [START_REF] Zeng | A comparative study on sufficient conditions for Takagi-Sugeno fuzzy systems as universal approximators[END_REF]. But it was clearly shown that the number of fuzzy rules increases as the approximation error tends to zero. So it becomes difficult to make use of the universal approximation property of T-S fuzzy modelling for practical purposes. Moreover, if the number of rules is bounded, the resulting set of functions is nowhere dense in the space of approximated functions [START_REF] Tikk | Polytopic and TS Models are Nowere Dense in the Approximation Model Space[END_REF]. These conflicting objectives have motivated researchers to find a balance between the specified accuracy and the computational complexity of the resulting fuzzy model.

In [START_REF] Al-Hadithi | New Methods for the Estimation of T-S Model Based Extended Kalman Filter and its Applications to Optimal Control for Nonlinear Systems[END_REF] and [START_REF] Jiménez | Improvement of Takagi-Sugeno Fuzzy Model for the Estimation of Nonlinear Functions[END_REF] new and efficient approaches are presented to improve the local and global estimation of T-S model. The first approach uses the minimum norm method to search for an exact optimum solution at the expense of increasing complexity and computational cost. The second one is a simple and less computational method, based on weighting of parameters.

The rest of the paper is organized as follows. A description of feedback linearization methodology is presented in section 2. In section 3, the estimation of T-S fuzzy model is discussed. In section 4, the design of an optimal LQR controller is developed. In section 5, an example of an inverted pendulum is illustrated to demonstrate the validity of the proposed approach. The conclusions of the validity of the proposed estimation approach and the robustness of the LQR are explained in section 6.

Feedback Linearization

Feedback linearization is an approach to nonlinear control design which has attracted a great deal of research interest. The central ideal is to algebraically transform a nonlinear system dynamics into (fully or partly) linear one so that linear control techniqus can be applied. This differes entirely from conventional linearization (i.e. Jacobian linearization) in that feedback linearization is achieved by exact state transformation and feedback, rather than by linear approximations of the dynamics. Supposing that the system can be modeled as follows:

y (n) (t) = f n (y(t), ẏ(t), • • • , y (n-1) (t)) + b n (y(t), ẏ(t), • • • , y (n-1) (t))u(t) (1)
Let us define the following vector:

x = y(t) ẏ(t) • • • y (n-1) (t)
T ∈ ℜ n

The system can be modeled as follows:

ẋn (t) = f n (x(t)) + b n (x(t))u(t) ( 2 )
For systems which can be expressed in the controllable and canonical form, the control action can be described as:

u(x) = 1 b n (x) (s(x) -f n (x)) (3)
Substituting (3) in (2) to obtain the feedback system:

ẋn (t) = f n (x) + b n (x)u(t) = f n (x) + b n (x) 1 b n (x) (s(x) -f n (x)) = s(x) (4)
Let us suppose the following control law:

s(x) = Sx = [-s 0 -s 1 • • • -s n-2 -s n-1 ]x (5)
The feedback system becomes a linear one with a characteristic polynomial p(λ) where the s i are chosen so that the the characteristic polynomial:

p(λ) = λ n + s n-1 λ n-1 + • • • + s 1 λ + s 0
The first inconvenient is the need to know explicitly the nonlinear system model, which in general can not be obtained because of the difficulty of identifying their parameters. Secondly, an adequate choice of the characteristic polynomial should be made. A possible solution for the first problem is by using the T-S fuzzy model for estimating nonlinear systems.

S (i1...in) : If z 1 is M i1 1 and . . . z f is M in n then ẋ = a (i1...in) 0 + A (i1...in) n x + B (i1...in) n u where A (i1••••••in) n = [a n1 (x) a n2 (x) • • • a nn (x)]
Applying this model, the identification of its parameters can be realized using the proposed estimation methods which will be described in section 3. We suppose that we have a first estimation of the T-S model parameters. In order to obtain such estimation, the classical least square method can be used around the equilibrium point. The objective is to obtain a local approximation of the system.

ẋn = a 0 o + a 0 n1 x 1 (x) + a 0 n2 x 2 (x) + • • • + a 0 nn x n (x) + b 0 n u (6)
A design of an optimal controller is carried out to show the effectiveness of the proposed estimation methods. The well known LQR controller is applied to solve the second problem proposed in this work which is how to choose the optimal characteristic polynomial p(λ) (see section 4).

Estimation of Fuzzy T-S Model's Parameters

An interesting method of identification is presented in [START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF]. The idea is based on estimating the nonlinear system parameters minimizing a quadratic performance index. It is based on the identification of functions of the following form:

f : ℜ n -→ ℜ y = f (v 1 , v 1 , . . . , v n )
Each IF-THEN rule for an n th order system can be written as follows:

S (i1...in) : If v 1 is M i1 1 and . . . v n is M in n then (7) ŷ = p (i1...in) 0 + p (i1...in) 1 v 1 + p (i1...in) 2 v 2 + . . . + p (i1...in) n v n
where the fuzzy estimation of the output is:

ŷ = r1 i1=1 . . . rn in=1 w (i1...in) (v) p (i1...in) 0 + p (i1...in) 1 v 1 + . . . + p (i1...in) n v n r1 i1=1 . . . rn in=1 w (i1...in) (v) (8)
where,

w (i1...in) (v) = n l=1 µ li l (v l )
being µ jij (v j ) the membership function that corresponds to the fuzzy set M ij j . Let {v 1k , v 2k , . . . , v nk , y k } be a set of input/output system samples. The parameters of the fuzzy system can be calculated as a result of minimizing a quadratic performance index:

J = m k=1 (y k -ŷk ) 2 = ∥Y -XP ∥
X = β (1...1) 1 β (1...1) 1 v 11 . . . β (1...1) 1 v n1 . . . β (r1...rn) 1 . . . β (r1...rn) 1 v n1 β 1...1 m β (1...1) m v 1m . . . β (1...1) m v nm . . . β (r1...rn) m . . . β (r1...rn) m v nm and β (i1...in) k = w (i1...in) (v k ) r1 i1=1 . . . rn in=1 w (i1...in) (v k ) (11)
If X is a matrix of full rank, the solution is obtained as follows:

J = ∥Y -XP ∥ 2 = (Y -XP ) T (Y -XP ) ∇J = X T (Y -XP ) = X T Y -X T XP = 0 (12) P = (X T X) -1 X T Y
In the case when the matrix X is not of full rank, an effective approach with few computational effort, based on the well known parameters' weighting method [START_REF] Al-Hadithi | Variable Structure Control with Chattering Reduction of a Generalized T-S Model[END_REF], [START_REF] Al-Hadithi | New Methods for the Estimation of T-S Model Based Extended Kalman Filter and its Applications to Optimal Control for Nonlinear Systems[END_REF] and [START_REF] Jiménez | Improvement of Takagi-Sugeno Fuzzy Model for the Estimation of Nonlinear Functions[END_REF]. This method can also be used for parameters tuning of T-S model from local parameters obtained through the identification of a system in an operating region or from any physical input/output data. We suppose that we have a first estimation

P 0 = [p 0 0 p 0 1 p 0 2 . . . p 0 n ] T
of the T-S model parameters. In order to obtain such an estimation, the classical least square method can be used around the equilibrium point. This first approximation can be utilized as reference parameters for all the subsystems. Then, the parameters' vector of the fuzzy model can be obtained minimizing:

J = m k=1 (y k -ŷk ) 2 + γ 2 r1 i1=1 . . . rn in=1 n j=0 (p 0 j -p (i1...in) j ) 2 = ∥Y -Xp∥ 2 + γ 2 ∥P 0 -p∥ 2 = Y γP 0 - V γI p 2 = ∥Y a -V a p∥ 2 (13)
where p 0 = [P 0 P 0 . . . P 0 ] T r1.r2...rn

In this case, the factor γ represents the degree of confidence of the parameters initially estimated [START_REF] Al-Hadithi | Variable Structure Control with Chattering Reduction of a Generalized T-S Model[END_REF], [START_REF] Al-Hadithi | Variable Structure Control with chattering elimination and guaranteed stability for a generalized T-S model[END_REF], [START_REF] Al-Hadithi | A new approach to fuzzy estimation of Takagi-Sugeno model and its applications to optimal control for nonlinear systems[END_REF] and [START_REF] Jiménez | Improvement of Takagi-Sugeno Fuzzy Model for the Estimation of Nonlinear Functions[END_REF].

Design of an Optimal Controller

A design of an optimal controller based LQR controller is applied to show the validity of the proposed estimation methods. The characteristic polynomial of the feedback system is calculated for this system so that finally, the parameters of this polynomial is used for the definition of s(x). In order to calculate the control coefficients, any control design methodology through state feedback control can be applied.

Optimal selection of closed loop poles will lead to a trade-off between speed of dynamic response and control effort. Together with the proposed estimation method, LQR might be an appropriate choice. With this regulator we guarantee the stability and gaining a balance between static and dynamic behavior of the system with admissible control actions. Selection of values of Q and R matrices determines the dynamic speed of the controller as well as amplitudes of state variables and control signals.

ẋ = Ax + Bu x ∈ ℜ n , u ∈ ℜ m , A ∈ ℜ n×n , B ∈ ℜ n×m
The objective is to find the control action u(t) to transfer the system from any initial state x(t 0 ) to some final state x(∞) = 0 in an infinite time interval, minimizing a quadratic performance index of the form:

J = ∞ t0 (x T Qx + u T Ru)dt
where Q ∈ ℜ n×n is a symmetric matrix, at least positive semidefinite one and R ∈ ℜ m×m is also a symetric positive definite matrix. The optimal control law is then computed as follows:

u(t) = -Kx(t) ( 1 4
)

K = -R -1 B T L (15)
where the matrix L ∈ ℜ n×n is a solution of the Riccati equation:

0 = -Q + LBR -1 B T L -LA -A T L
The feedback system becomes:

ẋ = Ax + Bu = Ax -BKu = (A -BK)x
and the characteristic polynomial is:

p(λ) = det(λI -A + BK)
It is aimed in this work to design an optimal LQR controller for the local approximation in [START_REF] Hseng | Fuzzy Bilinear Model and Fuzzy Controller Design for a Class of Nonlinear Systems[END_REF]. In order to obtain such estimation, the classical least square method can be used around the equilibrium point.

Application of the Proposed Control Algorithm for T-S Systems

The new proposed control algorithm can be considered as a two level one: the first level includes the calculation of the resulting matrices a 0 (x), A(x) and B(x) from the above identified fuzzy system and the second one is obtaining the control action mentioned above in [START_REF] Al-Hadithi | Variable Structure Control with chattering elimination and guaranteed stability for a generalized T-S model[END_REF]. The control law becomes:

u(x) = 1 b n (x) (SX -a 0n (x) -A n (x)x) = - 1 b n (x) (a 0n (x) + (A n (x) -V )x) (16)

Illustrative Example

Consider the problem of estimating an inverted pendulum using the above mentioned estimation methods. The inverted pendulum can be represented as follows:

θ = gsinθ -cosθ( u+ml θ2 sinθ M+m ) l( 4 3 -mcos 2 θ M+m ) (17)
Where θ denotes the angular position (in radians) deviated from the equilibrium position (vertical axis) of the pendulum and θ is the angular velocity, g(gravity acceleration), M(mass) of the cart=1 kg, m(mass) of the pole=0.1 kg, l is the distance from the center of the mass (m) of the pole to the cart=0.5 m. Assuming that x 1 = θ and x 2 = θ, then (17) can be rewritten in state space form as follows:

x 1 = θ ẋ1 = x 2 x 2 = θ ẋ2 = gsin(x 1 ) -cos(x 1 )( u+mlx 2 2 sin(x1) M+m ) l( 4 3 -mcos 2 (x1) M+m ) (18)
Firstly, the model of the inverted pendulum is estimated in three operation points for both the angle and its derivative. The universe of discourse of the angle is [ -π 4 , π 4 ] rad. and the one of the angular velocity is [-5, 5] rad seg . Both MFs for the angle x and its derivative ẋ are shown in figures 1 y 2 respectively. If we apply the T-S method directly to this example, then the condition number of the matrix X is 3.4148e +015, which shows clearly a non reliable result. Using the parameters' weighting method with weighting factor γ = 0.01, the inverted pendulum fuzzy model can be represented as follows:

S 11 : If x 1 is M 1 1 and x 2 is M 1 2 then ẋ2 = -8.2354 + 3.3866x 1 -0.2040x 2 -1.0443u S 12 : If x 1 is M 1 1 and x 2 is M 2 2 then ẋ2 = -8.3313 + 3.0983x 1 + 0.0000x 2 -1.0443u
. . .

S 33 : If x 1 is M 3 1 and x 2 is M 3 2 then ẋ2 = 8.0784 + 3.8914x 1 -0.2297x 2 -1.0903 (19)
By using the identification with the classical minimum square method in an interval close to the equlibrium point

x 1 ∈ [ -π 4 , π 4 ] x 2 ∈ [-2.5, 2.5]
The linear model of the system in this interval is:

ẋ2 = 0.0092 + 15.2665x 1 -0.0000x 2 -1.4187u
As it has been mentioned before, the inestable equilibrium point which will be also the objective of control is x 1 = x 2 = u = 0. In order that this point becomes an equlibrium one in the fuzzy model, a 22 0 should be zero. The resultant mean square error from this approximation is 0.0013. In this case, the condition number of the extended matrix X a becomes 1.4569e +004, thus improving the reliability of results.

As it can be observed, the results obtained through the parameters weighting method are always better than with the original T-S method. In fact, the proposed method correspond to the best possible result using TS method when the interval covering tends to totality, but this limit is not achievable since the problem would not be solvable any more.

Consider the problem of stabilizing and balancing of swing up of an inverted pendulum. The control of this system is a widely used performance measure of a controller, since this system is unstable and highly nonlinear. The objective is to maintain the inverted pendulum upright with θ despite small disturbances due to wind or system noises. The proposed optimum LQR is applied to calculate the desired parameters of the characteristic polynomial p(λ) defined in the previous rules, ẋ1 ẋ2 = 0 1 15.2665 0

x 1 x 2 + 0 -1.4187 u
minimizing the following performance index:

J = ∞ t0 (100x 2 1 + 10x 2 2 + u 2 )dt
and the computed control becomes u(t) = -Kx(t) = 25.4509 6.7734 x(t)

The matrix of the feedback system becomes:

A R = A -BK = 0 1 -20.8408 -9.6095
and its characteristic polynomial is:

p(λ) = λ 2 + 9.6095λ + 20.8408
and the eigen values are:

λ 1 = -6.3029 λ 2 = -3.3065
Thus, the control becomes:

u(x) = 1 b n (x) (a 0n (x) + (A n (x) -S)x) ( 2 0 )
where: S = -20.8408 -9.6095

Figure 3 shows the transient response of the inverted pendulum controlled by the proposed LQR subjected to an initial condition. The results obtained show that the system is stabilized by applying the proposed LQR. Moreover, we have proven that the LQR is robust and invariant against measurement noise. If we suppose that the angle sensor introduces a noise of σ = 1 • , it can be observed that the output is only affected by σ = 0.04 • as shown in figure 4. Figure 5 shows several trajectories in state space form of the system for several initial conditions. A feedback linearization is used to nonlinear control design for fuzzy systems described by T-S. The main objective is to obtain an improved performance of highly non-linear unstable systems. In this work, the well known weighting parameters approach is proposed to optimize local and global approximation and modelling capability of T-S fuzzy model to improve the choice of the performance index and minimize it. An inverted pendulum mounted on a cart is chosen to evaluate the robustness, effectiveness, accuracy and remarkable performance of proposed estimation approach. The results show that the proposed approach is less conservative than those based on (standard) T-S model. Simulation results indicate the potential, simplicity and generality of the estimation method and the robustness of the LQR algorithm. We prove that the proposed estimation algorithm converge very fast, thereby making them very practical to use.

2 ( 9

 29 ) where Y = y 1 y 2 . . . y m T P = p (1...1) 0 p (1...1) 1 p (1...1) 2 . . . p (1...1) n . . . p (r1...rn) 0 . . . p (r1...rn) n T[START_REF] Tikk | Polytopic and TS Models are Nowere Dense in the Approximation Model Space[END_REF]

Fig. 1 :Fig. 2 :

 12 Fig. 1: Membership functions for the angle x of the inverted pendulum Fig. 2: Membership functions for the angular velocity of the inverted pendulum

Fig. 3 :Fig. 4 :Fig. 5 :

 345 Fig. 3: Transient response of the inverted pendulum by the fuzzy LQR Fig. 4: Transient response of the inverted pendulum by the fuzzy LQR subjected to measurement noise

Acknowledgements

This work is funded by Spanish Ministry of innovation and Science (ARABOT project DPI2010-21247-C02-01).


Eleni Vrochidou

Petros Alvanitopoulos

Ioannis Andreadis

Anaxagoras Elenas

Katerina Mallousi

HHT-Based Artificial Seismic Accelerograms Generation

Keywords: Hilbert-Huang Transform, artificial spectrum-compatible seismic accelerograms, frequency components

A new efficient approach for generating spectrum-compatible seismic accelerograms is proposed. It is based on the Hilbert-Huang Transform (HHT); one natural seismic accelerogram is decomposed into frequency and amplitude components. The components are appropriately modified to synthesize the artificial seismic accelerogram that appears to have compatible acceleration spectrum with the natural seismic accelerogram. The HHT is an adaptive signal processing technique for analyzing nonlinear and non-stationary data such as seismic accelerograms. With HHT a seismic accelerogram is decomposed into a finite and small set of components. These components have well defined instantaneous frequencies, estimated by the first derivative of the phase of the analytic signal. The method is tested using twenty natural seismic records and a comparison with two established methodologies is provided.

Introduction

Earthquake engineering is the scientific field concerned with the study of the behavior of structures subject to seismic loading. Earthquake accelerograms (seismic signals) are required in order to simulate the response of structures. Disaster scenarios need an extensive set of seismic signals and although a large database of recordings of seismic excitations exists, for many regions there is lack of actual acceleration time-history records [START_REF] Carballo | Probabilistic seismic demand analysis: spectrum matching and design[END_REF]. Artificial spectrum-compatible accelerograms are widely used in the dynamic analysis of structures. The proposed method is based on HHT for generating spectrum-compatible earthquake accelerograms. It is well-known that the HHT is an appropriate signal processing technique for analyzing non-stationary and nonlinear signals such as seismic signals [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Alvanitopoulos | Interdependence between seismic parameters and damage indicies based on Hilbert-Huang Transform[END_REF].

The HHT decomposes the earthquake accelerogram into a finite number of Intrinsic Mode Functions (IMFs), and provides an energy-frequency-time distribution, the Hilbert spectrum. The extraction of the IMFs is based on the local characteristics of the seismic signal and not on a priori assumptions. Unlike traditional signal processing techniques, in HHT the frequency is defined through differentiation and there is no uncertainty principle limitation on time or frequency resolution from the convolution pairs based on a priori bases. Thus, HHT provides sharper energy-frequency-time distribution [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. Moreover, the finite number of the extracted IMFs is adequate to describe any seismic signal and therefore reduces the computational burden of the algorithm.

Previously reported methods in the literature employ HHT to generate spectrumcompatible accelerograms [START_REF] Ni | Application of Hilbert-Huang transform in generating spectrum-compatible earthquake time histories[END_REF][START_REF] Vrochidou | Correlation between seismic intensity parameters of HHT-based synthetic seismic accelerograms and damage indices of buildings[END_REF]. In method [START_REF] Ni | Application of Hilbert-Huang transform in generating spectrum-compatible earthquake time histories[END_REF] six natural records are required to describe the high and low frequency areas of the target design spectrum. The artificial spectrum-compatible seismic signal is obtained by solving an optimization problem. The difference between the target spectrum and the response spectrum of the artificially generated seismic signal is treated as the cost function. In contrast to the aforementioned method, in the proposed technique only one initial earthquake signal is required. Moreover, it can be of any intensity and thus no limitations are imposed. Another method is proposed in [START_REF] Vrochidou | Correlation between seismic intensity parameters of HHT-based synthetic seismic accelerograms and damage indices of buildings[END_REF].Artificial seismic signals are generated and submitted to a correlation study. The correlation study demonstrates that the artificially generated seismic signals share similar seismic parameter values to those of natural accelerograms.

The proposed methodology is applied to twenty natural seismic accelerograms. In order to reveal the effectiveness of the proposed method, a comparison with two wellestablished algorithms for generating spectrum-compatible accelerograms is presented. The first is Gasparini and Vanmarcke algorithm [START_REF] Gasparini | SIMQKE -A program for artificial motion generation, user's manual and documentation[END_REF] implemented by SIMQKE software [START_REF] Gelfi | SIMQKE_GR version 2.7[END_REF]. The second is Hallodorson and Papageorgiou algorithm [START_REF] Halldorsson | Calibration of the specific barrier model to earthquake of different tectonic regions[END_REF] implemented by SeismoArtif artificial earthquake accelerograms generator [START_REF]Computer program SeismoArtif, version 1.0.0[END_REF]. Results demonstrate that the proposed method is as reliable as the methodologies implemented by the two established Computer-Aided Design (CAD) tools.

2

The Hilbert-Huang Transform

The HHT was introduced by Huang et al [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] and comprises two parts: the Empirical Mode Decomposition (EMD) and the Hilbert Spectral Analysis (HSA).

Empirical Mode Decomposition (EMD)

The EMD can decompose complex data into a finite and small number of IMFs. The method is considered to be valid as it is an adaptive method and since the decomposition is based on the local characteristics of the signal, it is appropriate for nonstationary and non-linear signals [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Alvanitopoulos | Interdependence between seismic parameters and damage indicies based on Hilbert-Huang Transform[END_REF]. The decomposition process for a seismic signal is described in Figure 1 [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Alvanitopoulos | Interdependence between seismic parameters and damage indicies based on Hilbert-Huang Transform[END_REF].

Hilbert Spectral Analysis (HSA)

The HSA applies the Hilbert Transform to each IMF so as to obtain the instantaneous frequency data and construct the energy-frequency-time distribution. For all components, the Hilbert Transform is:

(t) c j (t) y j dτ τ t (τ) c Ρ j π 1 (t) j y ∫ +∞ ∞ - - = (1)
where denotes the Cauchy Principal Value. The IMF components are defined by the equation:

Ρ (2) (t) j (t)cosθ j a ) (t) j iθ (t)e j Re(a (t) j c = =
where is the real part. The initial signal is then written as: Re()

(3) ∑ ∫ = = n 1 j j j (t)dt)) 2πf (t)cos( a Re( X(t)
The residue is not included to the above equation because it is either a monotonic function or a constant. Equation [START_REF] Alvanitopoulos | Interdependence between seismic parameters and damage indicies based on Hilbert-Huang Transform[END_REF] shows that the amplitude and frequency can be expressed as functions of time and can strengthen the time-frequency distribution of the amplitude. This distribution is called the Hilbert spectrum.

(t) r n t) H(ω,

3

Response Spectrum and Strong Motion Duration

Response Spectrum

The response spectrum is an essential tool in earthquake engineering and is the plot of the peak response acceleration, displacement or velocity of a single-degree-of-freedom oscillator in dependence of its fundamental period under the same seismic excitation. Figure 2 demonstrates the concept of the earthquake response spectrum.

Strong Motion Duration (SMD)

A major issue of the spline fitting arises at the boundaries since cubic splines have wide swings if they are left unattended [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. To eliminate end effects, only the strong motion duration of the seismic excitation is considered in the proposed methodology.

The strong motion duration (SMD) is the time of the earthquake where the most seismic energy is released. The SMD after Trifunac/Brady [START_REF] Trifunac | A study on the duration of strong earthquake ground motion[END_REF] is defined as time elapsed between 5% and 95% of the Husid diagram [START_REF] Husid | Analisis de terremoros: analisis general[END_REF] and is defined as follows:

0.05 0.95 0.90

T T T - = (4)
where T 0.90 is the SMD and T 0.95 and T 0.05 are the time elapsed at the 95% and 5% of the Husid diagram, respectively.

Overview of the Proposed Method

According to the proposed methodology, the original signal is subjected to the HHT so as to obtain the n amplitude and frequency components. The frequency components of the SMD are submitted to a smoothing procedure and addition of random Gaussian noise. The amplitude components of the SMD are also modified. The artificially generated accelerogram is derived by reconstructing the seismic signal from the new modified instantaneous frequency and amplitude components:

t) AGA(x, jNew ω jNew a (5) ∑ = ∫ = n 1 j (t)dt ω j jΝew jΝew (t))e (a Re t) AGA(x,
Each new amplitude component is given by:

, (t) a x (t) a jΟld k jΝew = n ,..., 2 1, k = (6)
where, is the number of extracted IMFs of the initial seismic signal, x is a vector containing scaling parameters of values randomly distributed in the interval [0.8-1.3] and is the initial amplitude value at time t. The new frequency content is obtained after a two step procedure performed to the frequency components of the natural seismic signal extracted from the EMD. This procedure is applied only to the IMFs that have in the SMD, mean frequency value lower than 5 Hz. This usually occurs after the third IMF and above. The rest of the IMFs remain unaffected. The response spectrum acceleration is determined in period range 0.2 s to 4 s (0.25 Hz to 5 Hz). This range covers the majority of low-, mid-and high-buildings [START_REF] Bommer | Displacement spectra for seismic design[END_REF]. The first three IMFs are of high frequency, thus of low eigenperiod, so there is no meaning of changing them as they do not contribute to the approximation of the spectrum at the eigenperiods of our interest.

The first step of the procedure employs a smoothing process to the SMD of the appropriate frequency IMFs by utilizing a moving average filter. At a second stage, Gaussian noise is added to create new frequency components that could resemble to the frequency components that were derived from the initial signal. Figure 3 shows the block diagram of the proposed method.

Figures 4(a For the specific earthquake signal the total duration is 24.44 s while the strong motion duration is between 4.15 s and 12.67 s.

Empirical Results

The experimental results presented in this section include a comparison of the proposed method with similar state of the art techniques that are widely used for generating artificial spectrum-compatible accelerograms. Database, numerical results and system setup used in the validation process are described in the following sections.

Experimental Setup

In order to assess the performance of the proposed method, an experimental comparison is performed with respect to two other supervised approaches, invoked by the programs SIMQKE [START_REF] Gelfi | SIMQKE_GR version 2.7[END_REF] and SeismoArtif [START_REF]Computer program SeismoArtif, version 1.0.0[END_REF]. SIMQKE is based on Gasparini and Vanmarche algorithm [START_REF] Gasparini | SIMQKE -A program for artificial motion generation, user's manual and documentation[END_REF]. It estimates a power spectral density function from a specific response spectrum and generates the amplitudes of sinusoidal signals that have phase angles randomly distributed between 0 and 2π. The SeismoArtif, provides the Synthetic Accelerogram Generation and Adjustment method, based on Hallodorson and Papageorgiou algorithm [START_REF] Halldorsson | Calibration of the specific barrier model to earthquake of different tectonic regions[END_REF]. From a synthetic accelerogram simulated by the user, the spectrum-compatible artificial accelerogram is defined by adapting its frequency content using the Fourier Transform. The Root Mean Square Error (RMSE) between the response spectrum of the generated accelerogram and the spectrum of the natural accelerogram is calculated in all experiments:

∑ = - = n 1 j 2 target generated ) SA SA ( Ν 1 RMSE (7)
where, N is the number of samples. The calculation of the RMSE is utilized for the numerical validation of the compatibility between the response spectra of the original and generated seismic accelerograms.

Performance Evaluation

Table 1 includes all information regarding the seismic events [START_REF]PEER Ground Motion Database[END_REF], that are selected to test the method, and the values of RMSE for the proposed method (PM), the SIMQKE (M1) and the SeismoArtif (M2). It can be observed that all RMSE values are very close in all methods. However, the PM leads to spectrum compatible artificial seismic accelerograms with lower RMSE in most experiments. The PM can be directly compared with the two wellknown procedures and be advantageous since is not an iterative method. The RMSE calculated for the two standard methods is the average of ten measurements derived from different number of iteration [START_REF] Carballo | Probabilistic seismic demand analysis: spectrum matching and design[END_REF][START_REF] Alvanitopoulos | Interdependence between seismic parameters and damage indicies based on Hilbert-Huang Transform[END_REF][START_REF] Vrochidou | Correlation between seismic intensity parameters of HHT-based synthetic seismic accelerograms and damage indices of buildings[END_REF][START_REF] Gelfi | SIMQKE_GR version 2.7[END_REF][START_REF]Computer program SeismoArtif, version 1.0.0[END_REF][START_REF] Husid | Analisis de terremoros: analisis general[END_REF][START_REF]PEER Ground Motion Database[END_REF][START_REF]EC8, Eurocode 8: Design of Structures for Earthquake Resistance-Part 1: General Rules, Seismic Actions, and Rules for Buildings[END_REF], 20 and 25) while the RMSE for the PM is the average of ten measurements of only one iteration each. In Table 2 it is recorded the RMSE for six experiments with the PM, and the number of iterations that need both M1 and M2 in order to obtain approximately the same range of values for the RMSE. The PM performs one iteration to achieve that range of RMSE, while at least three iterations are required for M1 and over five for M2. Figure 5 shows the natural seismic signal of Victoria, and the generated seismic signal with the PM, M1 and M2. The artificial signals generated by the two programs and the PM are of the same duration and compatible with the same target spectrum. SIMQKE generates artificial signals by varying the random phase angles of the sinusoidal signals and all generated signals are restricted to have a uniform shape that is provided by a predetermined envelop. In SeismoArtif, the generation of the synthetic accelerogram starts from a Gaussian white noise that is multiplied by the envelop shape suggested by Saragoni and Hart [START_REF] Saragoni | Simulation of artificial earthquakes[END_REF]. Both methods derive accelerograms that follow predetermined envelop shapes, thus, accelerograms generated with either of the two methods tend to resemble among them.

In order to demonstrate the ability of the PM to generate multiple different accelerograms, the standard deviation (SD) is used as a metric. Figure 6 demonstrates the SD between ten artificially generated seismic signals with the PM and the mean of them, starting from the natural record of Victoria, in contrast to the SD of ten artificially generated seismic signals with M1 and M2 and the mean of them. During the SMD of the signal, the SD is greater with the PM rather than with M1 and M2. The mean spectrum is the average of ten artificially generated signals. The mean spectrum of the PM is evaluated from ten measurements with different amount of smoothing by varying the span of the moving average filter and adding more noise to the components. The mean spectrum of M1 and M2 stem from ten measurements with different number of iterations.

Finally, all generated artificial seismic accelerograms with the PM follow the restrictions of Eurocode 8 [START_REF]EC8, Eurocode 8: Design of Structures for Earthquake Resistance-Part 1: General Rules, Seismic Actions, and Rules for Buildings[END_REF]. According to Eurocode 8, in the range of periods between 0.2T 1 and 2T 1 (T 1 is the fundamental period of the structure) no value of the mean 5% damping elastic spectrum, calculated from all time histories, should be less than 90% of the corresponding value of the 5% damping elastic response spectrum. In Figure 7(b) are demonstrated the 90% of the target response spectrum and the mean spectra of four artificially generated seismic signals with the PM.

Conclusions

Artificial spectrum-compatible accelerograms are used instead of natural acceleration records for the dynamic analysis of many critical structures for sites with no strong ground motion data. A new method for generating spectrum-compatible seismic accelerograms based on the HHT is proposed. One randomly selected seismic event is decomposed in amplitude and frequency components. The amplitude components of the SMD are scaled and the appropriate frequency components of the SMD are subjected to smoothing and addition of Gaussian noise. The PM is flexible, since it requires only one seismic event, that can be from any region globally and of any intensity. This option makes the algorithm faster without affecting its accuracy by choosing seismic accelerograms with a definite number of points. Moreover, in order to analyze the non-stationary data only a small number of IMFs is required.

Comparison with two other state of art methods, one based on Gasparini and Vanmarcke algorithm (calculated by SIMQKE) and the other based on Hallodorson and Papageorgiou algorithm (calculated by SeismoArtif), has shown the satisfactory performance of the PM. The introduction of SD revealed that the PM can derive multiple different accelerograms. The proposed method is not an iterative method and all generated signals follow the restrictions of the Eurocode 8.

Fig. 1 .

 1 Fig. 1. Flow chart of the EMD.

Fig. 2 .

 2 Fig. 2. Concept of the earthquake response spectrum.

Fig. 3 .

 3 Fig. 3. Block diagram of the proposed method.

Fig. 4 .

 4 Fig. 4. (a) The 3rd and (b) 4th IMF of Victoria, Mexico time history record after HSA and (c), (d) after smoothing and addition of Gaussian noise.

  ) and 4(b) show two frequency components of the earthquake in Victoria Mexico. The smoothing and addition of Gaussian noise are performed to the SMD of the earthquake signal and the new components are generated and depicted in Figures 4(c) and 4(d).

Fig. 5 .

 5 Fig. 5. (a) Victoria, Mexico time history record (b) artificially generated accelerogram with the PM (c) artificial seismic signal derived from M1 (d) artificial seismic signal derived from M2.

Figure 7 (

 7 Figure7(a) shows the target spectrum and the mean spectra of the proposed and the two established methods evaluated with the programs SIMQKE and SeismoArtif. All spectra are compatible with the target spectrum and have slight differences (Table1).

Fig. 6 .

 6 Fig. 6. Standard deviation between ten artificially generated signals (compatible with the Victoria spectrum) and the mean of them, using (a) the PM and SIMQKE (M1) (b) the PM and SeismoArtif (M2).

Fig. 7 .

 7 Fig. 7. (a) Target spectrum and mean spectra of the PM, M1 and M2 (b) target spectrum, spectra of four artificial generated signals with the PM and 90% of the target spectrum.

Table 1 .

 1 Numerical results of the method for 20 seismic events.

	Seismic event Country	Date Station / Component	RMSE PM	RMSE M1	RMSE M2
	Mammoth Lakes U.S.A. 05/27/80	FIS / 000	1.80•10 -2 1.46•10 -2 5.09•10 -2
	Dinar	Turkey 01/10/95	Dinar / 090	4.94•10 -2 2.45•10 -2 4.89•10 -2
	Victoria	Mexico 06/09/80 Station 6604 / 045	2.04•10 -2 3.08•10 -2 3.43•10 -2
	Gazli	Uzbekistan 05/17/76	Karakyr / 000	3.79•10 -2 1.38•10 -2 3.60•10 -2
	Tabas	Iran	09/16/78	Tabas / LN	9.12•10 -2 1.17•10 -1 7.81•10 -2
	Spitak	Armenia 12/07/88	Gukasian / 000	1.06•10 -2 1.31•10 -2 4.69•10 -2
	Morgan Hill	U.S.A. 04/24/84 Station 1180 / 000	8.86•10 -3 3.84•10 -3 1.52•10 -2
	Patra	Greece 07/14/93 Station 00001 / NS 5.05•10 -2 1.31•10 -1 5.83•10 -2
	Pyrgos	Greece 03/26/93 Station 00001 / NS 4.07•10 -2 5.33•10 -2 7.28•10 -2
	Loma Prieta	U.S.A. 10/18/89	Santa Cruz / EW	2.60•10 -2 4.96•10 -2 3.46•10 -2
	Big Bear	U.S.A. 06/28/92 Lake-Civic Center / 270 2.94•10 -2 2.55•10 -2 5.29•10 -2
	Dinar	Turkey 01/10/95	Balikesir / 090	2.29•10 -4 3.10•10 -4 1.49•10 -2
	Erzican	Turkey 03/13/92	Erzican / NS	1.57•10 -2 2.66•10 -2 3.53•10 -2
	Friuli	Italy	05/06/76	Barcis / 270	1.43•10 -3 1.51•10 -3 1.10•10 -2
	Imperial Valley	U.S.A. 10/15/79 Bonds Corner / 230 5.68•10 -3 4.56•10 -3 3.25•10 -2
	Kobe	Japan 01/16/95	Fuk / 090	7.17•10 -3 3.31•10 -3 3.16•10 -2
	Loma Prieta	U.S.A. 10/18/89 Apeel 10 Skyline / 090 9.47•10 -3 2.52•10 -2 7.48•10 -2
	Mammoth Lakes U.S.A. 05/27/80 Convict Creek / 090 1.20•10 -2 9.88•10 -3 2.24•10 -2
	Northridge	U.S.A. 03/20/94 Anacapa Island / 270 1.24 •10 -3 3.97•10 -4 2.53•10 -2
	San Fransisco	U.S.A. 03/22/57 Golden Gate / 100	4.52•10 -3 8.96•10 -3 2.07•10 -2

Table 2 .

 2 RMSE toward number of iterations for the proposed method, SIMQKE and Seis-moArtif for six seismic events.

	Seismic event (Station)	Method	RMSE	Number of iterations
	Mammoth Lakes (FIS)	PM M1 M2	1.80•10 -2 5.76•10 -2 1.74•10 -2	1 11 19
	Dinar (Dinar)	PM M1 M2	4.94•10 -2 3.60•10 -2 4.58•10 -2	1 3 7
	Victoria (Station 6604)	PM M1 M2	2.04•10 -2 2.25•10 -2 2.77•10 -2	1 19 25
	Gazli (Karakyr)	P.M M1 M2	3.79•10 -2 4.12•10 -2 3.87•10 -2	1 3 11
	Tabas (Tabas)	PM M1 M2	9.12•10 -2 9.21•10 -2 9.27•10 -2	1 7 5
	Spitak (Gukasian)	PM M1 M2	1.06•10 -2 1.04•10 -2 1.08•10 -2	1 5 25

Table


Christos Makris
email: makri@ceid.upatras.gr

Michael Angelos Simos
email: asimos@ceid.upatras.gr

Novel Techniques for Text Annotation with Wikipedia Entities

Keywords: Text semantic annotation, Wikipedia entities, Semantic Data Linking on Web, Data Mining, Information Retrieval, Ontologies

Text annotation is the procedure of identifying the semantically dominant words of a text segment and attaching them with conceptual content information in their context. In this paper, we propose novel methods for automatic annotation of text fragments with entities of Wikipedia, the largest knowledge base online, a process commonly known as Wikification aiming at resolving the semantics of synonymous and polysemous terms accurately. The cornerstone of our contribution is a novel iterative Wikification approach, converging at optimal annotations while balancing high accuracy with performance. Our first two methods can be fine-tuned through a machine-learning technique over large homogenous data sets. Our experimental evaluation resulted in remarkable improvement over state-of-the-art Wikification approaches.

Introduction

Wikipedia is one of the largest online knowledge repositories, consisting of more than 4 million entities (also called lemmas or articles), built over more than with over 20 million contributing users. The online "crowdsourcing" nature of Wikipedia enables the creation and maintenance of knowledge entities in a quite diverse manner resulting in widespread and commonly accepted textual descriptions, through the consensus of a large number of people. Wikipedia inherits the main principles of web, elegantly combining textual content and cross-references with hyperlinks providing plenteous semantic information that can be exploited for representing the described entities.

Text annotation is the process of identifying the semantically dominant words (called spots or anchors) in a text segment and attaching them with additional information expressing conceptual content in their current context. Recent research works [START_REF] Cucerzan | Large-Scale Named Entity Disambiguation Based on Wikipedia Data[END_REF][START_REF] Ferragina | TAGME: on-the-fly annotation of short text fragments (by wikipedia entities)[END_REF][START_REF] Gabrilovich | Computing semantic relatedness using Wikipedia-based explicit semantic analysis[END_REF][START_REF] Han | Collective entity linking in web text: a graph-based method[END_REF][START_REF] Han | Named entity disambiguation by leveraging wikipedia semantic knowledge[END_REF][START_REF] Kulkarni | Collective annotation of Wikipedia entities in web text[END_REF][START_REF] Mihalcea | Wikify!: linking documents to encyclopedic knowledge[END_REF][START_REF] Milne | Learning to link with Wikipedia[END_REF][START_REF] Ratinov | Local and global algorithms for disambiguation to Wikipedia[END_REF] explored various approaches for efficient text annotation with Wikipedia entities. Text annotation is a fundamental preprocessing step of many information retrieval activities like semantic indexing, cross-referencing content on web, clustering, classification, for efficient rankings, text summarization, natural lan-guage generation, exploiting conceptual similarity of documents, and sentiment analysis.

Text annotation techniques face the entity disambiguation problem, where multiple candidate Wikipedia article annotations exist for a specific term, deriving by polysemy and synonymy, common properties of natural languages. Disambiguating to Wikipedia entities is similar with classic Word Sense Disambiguation tasks [START_REF] Navigli | Word Sense Disambiguation[END_REF], but also requires compatibility with additional information provided by the link structure of Wikipedia. The core of most existing methods involves the textual context and the collective agreement of other identified spots with every candidate disambiguation of a specific spot.

The aim of current work is the introduction of simple and accurate disambiguation methods for small texts, using Wikipedia as the underlying catalogue, thus focusing on the same kind of texts as in TAGME system [START_REF] Ferragina | TAGME: on-the-fly annotation of short text fragments (by wikipedia entities)[END_REF] and the approaches introduced in [START_REF] Meij | Adding semantics to microblog posts[END_REF][START_REF] Milne | Learning to link with Wikipedia[END_REF].

Literature Overview

Wikipedia disambiguation techniques can be distinguished into two main classes: local and global [START_REF] Ratinov | Local and global algorithms for disambiguation to Wikipedia[END_REF]. Local approaches exploit context information concerning only the specific spot, while global approaches rely on a collective agreement of all spots in a text segment and their disambiguation. Moreover, some approaches involve collection-level methods, performing a collective disambiguation of all entity references across a given document collection. Disambiguation based on Wikipedia entities has attracted many researchers [START_REF] Cucerzan | Large-Scale Named Entity Disambiguation Based on Wikipedia Data[END_REF][START_REF] Kulkarni | Collective annotation of Wikipedia entities in web text[END_REF][START_REF] Meij | Adding semantics to microblog posts[END_REF][START_REF] Mihalcea | Wikify!: linking documents to encyclopedic knowledge[END_REF][START_REF] Milne | Learning to link with Wikipedia[END_REF][START_REF] Ratinov | Local and global algorithms for disambiguation to Wikipedia[END_REF], with most prominent the approach proposed in [START_REF] Ferragina | TAGME: on-the-fly annotation of short text fragments (by wikipedia entities)[END_REF], where the authors focused on real time annotation of short texts existing on tweets, web search engine results etc, where textual context is quite limited. The TAGME system was extended in [START_REF] Makris | Improved Text Annotation with Wikipedia Entities[END_REF] by using information from knowledge resources such as WordNet [START_REF] Fellbaum | WordNet, an electronic lexical database[END_REF], and exploiting the PageRank values inherent in the Wikipedia pages. The system can be extended to handle other ontologies [START_REF] Hoffart | YAGO2: exploring and querying world knowledge in time, space, context, and many languages[END_REF] as well.

Proposed Methods

The touchstone of our proposed techniques is TAGME's voting scheme that assigns to each candidate sense of a spot a value, based on the collective agreement of this annotation from the other spots in the same text. Our first approach investigates some minor changes over TAGME's scheme. The second method proposes a more analytical model for the problem. Finally our iterative method, introduces a novel approach.

Spots Extraction

We use a simple yet very effective method for preprocessing the input text fragments for the extraction of candidate Wikipedia entities. Firstly, every text fragment is separated into tokens (words, punctuation, symbols, etc) to be grouped later, forming phrases when possible. The grouping is performed by matching as many adjacent tokens as possible with phrases (token groups) extracted from Wikipedia's links anchors. Complex rules have been applied allowing loose matching, handling minor details and further improving our results. The final step of spots extraction process is the removal of words not containing semantic information, (stopwords), symbols, punctuation etc. A more complex morphosyntactic analysis model is to be developed as future work.

Method 1

Let be an anchor to be annotated, and Pg( ) the list of candidate Wikipedia articles for this annotation, derived from occurrences of in Wikipedia hyperlinks as anchor to the above articles. A vote provided by a spot c to the annotation → with in Pg( ) is evaluated as:

( ) = ∑ ( , ) ( ) * ( | ) / | ( )|
where P( |c) is the prior probability of c is pointing to p l (also known as commonness). Commonness is pre-calculated by parsing the Wikipedia dataset and for each observed anchor-Wikipedia pair dividing the number of appearances that the anchor points to the corresponding article by the total appearances of the anchor. The quantity srel( , ) is a measure devised by us to depict the relatedness between two Wikipedia pages that is computed by taking into account the commonality of the incoming Wikipedia links to the two pages, as the cardinality of intersection divided by the maximum cardinality set, thus:

srel( , ) = | ( )∩ ( )| !"# ($ % &$,| ( )|)
where in(p) as in [START_REF] Fellbaum | WordNet, an electronic lexical database[END_REF][START_REF] Gabrilovich | Computing semantic relatedness using Wikipedia-based explicit semantic analysis[END_REF] denotes the set of Wikipedia pages pointing to page p.

The disambiguation weight of a probable annotation → is defined as:

" ' % & = ( " ) ( ) " ) *" '
The final annotation is selected based on maximum commonness score after a filtering step among Wikipedia articles within a voting score distance threshold τ from the maximum voting score.

We have tested relatedness formulas like [START_REF] Han | Collective entity linking in web text: a graph-based method[END_REF][START_REF] Milne | Learning to link with Wikipedia[END_REF], but the highest accuracy is achieved on this schema by our simple formula (srel), evaluating the incoming link intersection of two articles. A common observation throughout our experiments is that the heterogeneity of compared articles in size, results in imbalances of relatedness score, when using approaches like [START_REF] Han | Collective entity linking in web text: a graph-based method[END_REF] and [START_REF] Milne | Learning to link with Wikipedia[END_REF]. A better formula is to be developed as future work.

Method 2

Let Pg( ) be the set of candidate Wikipedia articles ( " ' j) annotations for each spot as anchor . The anchor-Wikipedia pairs that appear less than three times, or in less than 0.1% of total occurrences of the anchor are filtered. Let a text fragment contain k spots where i ϵ {1 , 2, …, k}.

The core logic of this method involves computing a global score for every probable combination of anchor senses represented as possible Wikipedia article annotations ( " ' ). There are

∏ % |,-" ' | . & / . = ∏ | | / .
probable combinations, that can be reduced by filtering more anchor-Wikipedia low probability pairs. Let Gscore be the global evaluation score of an annotation combination, as described above:

0 % " 1 , " 2 … " 4 & = ( ( 5 " ' , " 6 / 78 9.
/:.

8.

At first in "Method 2 (rel)" we evaluated a relatedness measure between two Wikipedia Articles:

, = 1 - %= >% $?@% &$, |?@ |&& - |?@ ∩ ?@ | A - =?@ |?@ |, |?@ |
proposed by Milne & Witten [START_REF] Milne | Learning to link with Wikipedia[END_REF] as the Bscore formula, where W denotes the size of Wikipedia. Still our simpler relatedness approach (srel) outperforms this commonly utilized relatedness measure yielding more accurate results, for most entities pairs having small incoming links intersection sets.

Finally, in "Method 2 comrel" we tested a more complex formula for Bscore:

5 % " ' , " 6 & = = % " ' , " 6 & = % " ' , " 6 & " ' | " 6 | 7
where P(p ai |a i ) is the Commonness of the anchor a i pointing the article p ai . This way BScore scales respectively with the commonality of the two articles combination. As a final step, again we perform filtering by keeping all candidate combinations achieving a Gscore up a threshold distance τ from the maximum Gscore. Among the filtering results, we select the combination maximizing the total commonness ∑ " ' , B C8D

. For very small threshold values (t), the selection of annotations yielding the highest Bscore is voted as the dominant. If the threshold is large, then a selection with the highest total commonness score is selected.

Threshold Optimization.

The most important fine tuning in our methods schemas involves the optimal selection of τ. We propose a method for optimizing τ threshold, given a Wikipedia entity manually annotated training set, having homogenous statistics with the to-beannotated dataset. This technique can also be used to determine experimentally the upper bound of performance for our methods (by training with our experimental dataset).

Based on this optimization technique we explore the threshold value intervals for which a selection with the optimal annotations passes through the final filtering step. Then we just find the values interval for τ, that belongs to the maximum number of each selection interval, for the selections that perform a correct annotation.

Thus, the problem of optimizing τ value can be deducted to a simple optimization problem of selecting a value that belongs to as more intervals as possible, by simply checking the values on intervals edges.

3.3

A Novel approach on text annotation.

Various model schemas have been introduced for the text annotation problem in previous works with most common the voting scheme and its many variations, to which we contributed during the previous part of this paper. One of the most important breakthroughs of the current work is the introduction of a novel method for semantic anchor annotation with Wikipedia entities, through an iterative approach, aiming at converging to an optimal candidate for each anchor as a result of constant improvement of the approximate solutions of each iteration. The principal intuition of this method, lies on the analysis of human semantic interpretation process of text segments containing polysemous terms. Such a process would involve an iterative evaluation procedure, until some termination criteria are met. For the most of the common text segments disambiguation, those criteria may be met before iterations are required, but in complex polysemous context, an evaluation of each candidate annotation may be necessary, consisting an iterative procedure until some criteria for a decision with a degree of certainty are met. The main logic difference of this method and method 2, which evaluates every possible annotation combination among anchors of a text segment, is the targeted evaluation of semantically meaningful combinations (as evaluated by some commonness and relatedness formulas), resulting a vast complexity reduction.

A candidate model of such a process may be epitomized by the following iterative method:

Let s 0, ,s 1 ,s 2 ,… s n be the spots of a text fragment, and P si 1 … P si m ϵ Pg(s i ) , the m candidate Wikipedia entity annotations of spot s i , where i ϵ {0,1,2,3…n}. Let Pg(s i ) the set of candidate Wikipedia entity annotations of spot s i . Each state is constituted of lists of candidate Wikipedia entity annotation of each spot sorted by a ranking criterion. At initial state, the candidate annotation lists are sorted by commonness in Wikipedia where p si 0, the most common entity of each spot s i ).

s So the disambiguation approximation during the initial step of the iteration is p s0 0 p s1 0 p s2 0 … p sn 0, equivalent with the commonness annotation approach. The iterative step involves sorting each of the spots list in descending order of its elements R metric values. Every iteration results in approximation improvement of the optimal annotation entities, by exploiting as initial seed the commonness of the candidate entities, and combining relatedness at each iteration step until convergence. The evaluated formula for the R metric is:

R( E ' x) = ∑ " ' , " 6 * " ' | * " 6 | 7

78D

Where: , =

|?@ F ∩?@ | = > G?@H F IG,$?@% &$

The iterative algorithm is terminated when one of the following convergence criteria is met:

─ A maximum number of iterations is reached. This is more of a security criterion, for handling the case of infinite iterations, due to other convergence criteria not being met. It is rarely activated, since in most cases, the rest of the convergence conditions are already met. The enforcement of large values results in more accurate results, relying more and more on relatedness for the cases where one of the other termination conditions is not yet satisfied. The application of smaller values, weights commonness more over relatedness, thus results in faster overall and worst case execution time, with an expense in accuracy. ─ There are zero disturbances during the last k iterations. This means that a convergence is due to happen, since no changes on the list rankings are observed. ─ The top m elements of each spots list maintain their position for k iterations. Thus the probability of disturbances is decreased. ─ Small difference of each of the spot lists elements values is observed between iteration steps, stabilizing the ranking and leading to convergence.

The convergence speed of the above criteria can be enhanced by filtering the lower k percentile elements of each lists spots, after m iterative steps, since our evaluation revealed a fairly decreased probability of the lowest k elements occupying top positions on the lists rankings, before some of the convergence criteria is satisfied. This behavior is intuitively explained by the fact that for the vast majority of cases, the ranking results tend to converge at optimal values through each iteration.

The real advantages of the third method, involve remarkable improvements compared with our two first contributions. Its iterative nature, allows balancing between speed and accuracy, with the appropriate fine tuning of the convergence conditions parameters. Flaccid convergence criteria lead to more mature results in terms of relatedness, approximating at most cases the accuracy of our second method. Tighter convergence criteria imply improved average and worst case time complexity, without significant sacrifices in accuracy. Finally we should note that convergence occurs within the first iterations in most cases with high accuracy, a fact that allows us conclude correlation between convergence speed ( in iterations) and the degree of certainty of the results accuracy, as is intuitively expected.

We exploited the above correlation among convergence speed and accuracy, by composing a post processing pruning step and evaluating the experimental results during our evaluation.

Experimental Evaluation

For the experimental evaluation of the methods we introduced above, we used the datasets of TAGME [START_REF] Ferragina | TAGME: on-the-fly annotation of short text fragments (by wikipedia entities)[END_REF], available online by the authors. Our principal focus involved the Wiki-Disamb30 data file consisting of 1.4M short texts randomly selected from Wikipedia pages, each one containing approximately 30 words and at least one ambiguous anchor. Common precision metrics were used as in [START_REF] Milne | Learning to link with Wikipedia[END_REF][START_REF] Kulkarni | Collective annotation of Wikipedia entities in web text[END_REF][START_REF] Fellbaum | WordNet, an electronic lexical database[END_REF], in order to evaluate in practice the performance of our algorithms.

In table 1 and figure 1, we present the precision of methods 1 and 2, after optimization of their parameters after training. The optimal result column of table 1, presents the maximum possible accuracy of the method, using optimal parameters, calculated through training with the entire evaluated dataset. Figure 2, shows the precision of method 3, by varying the number of max iterations parameter. Figure 3 presents the evaluation of method 3 pruning step (precision/recall of pruning). The evaluation of method 3 is presented separately from methods 1 and 2, since it doesn't involve a machine learning procedure.

Figure 1

We used an Ubuntu Server 12.04LTS on a quad core 2.9 GHz 64 bit PC utilizing 8GBytes of main memory. The proposed methods were implemented in Python 2.7.6 over the Wikipedia dump files instance enwiki-20130503, loaded on a PostgreSQL database.

We ran our methods on randomly selected subsets of Disamb30 dataset, since the time complexity of method 2 was beyond our time resources for a full scale exper ment, having in all cases similar results. The results, provided below are over a subset of more than 10K disambiguated spots. We also performed a simplified implement tion of TAGME disambiguation algorithm, as described in line, applying a similar threshold optimization technique of its Table 1 Table 1 and figure methods proposed, expressing their capacity to give the same senses/Wikipedia a ticles initially attached to the text fragments anchors. Throughout our experiments, we investigated the upper bound of our methods precision whi ble anchor of the dataset.

Our first method outperforms TAGME, yielding more accurate to 92.43%. This improvement is an outcome of the use of necessity of a better relatedness measur annotation algorithm schemas. The behavior of this method under our technique was also improved comparatively with TAGME, converging faster near optimal results.

The second method reached optimal of this method was by far larger than both our first method and TAGME. The rel and srel formulas were not as effective, yielding respectively 92. accuracy. But when using training set size for fast convergence to optimal results was remarkable. This fact was analyzed in depth, concluding that ranking quality pre-filtering list, with the correct that list in most cases. That leads to very small optimal threshold values ( thus fast convergence to the optimal ing commonness with relatedness by multi methods on randomly selected subsets of Disamb30 dataset, since the time complexity of method 2 was beyond our time resources for a full scale exper ment, having in all cases similar results. The results, provided below are over a subset isambiguated spots. We also performed a simplified implement tion of TAGME disambiguation algorithm, as described in [START_REF] Ferragina | TAGME: on-the-fly annotation of short text fragments (by wikipedia entities)[END_REF] as performance bas line, applying a similar threshold optimization technique of its ϵ value.

Figure 2

Table 1 and figure 1 present our detailed results on the precision of the expressing their capacity to give the same senses/Wikipedia a ticles initially attached to the text fragments anchors. Throughout our experiments, we investigated the upper bound of our methods precision while annotating every avail ble anchor of the dataset.

Our first method outperforms TAGME, yielding more accurate optimal . This improvement is an outcome of the use of srel formula, depicting the necessity of a better relatedness measure between two Wikipedia entities for similar annotation algorithm schemas. The behavior of this method under our τ-optimization technique was also improved comparatively with TAGME, converging faster near

The second method reached optimal accuracy of 94.61%, but the time complexity of this method was by far larger than both our first method and TAGME. The rel and srel formulas were not as effective, yielding respectively 92.07% and 92.98%, optimal accuracy. But when using comrel formula, the improvement in both accuracy and training set size for fast convergence to optimal results was remarkable. This fact was analyzed in depth, concluding that the analytic model of Method 2 creates a high filtering list, with the correct article of an annotation in the top of that list in most cases. That leads to very small optimal threshold values (τ) selection, thus fast convergence to the optimal τ. In conclusion, we note the necessity of weigh ing commonness with relatedness by multiplying as our comrel formula has proven methods on randomly selected subsets of Disamb30 dataset, since the time complexity of method 2 was beyond our time resources for a full scale experiment, having in all cases similar results. The results, provided below are over a subset isambiguated spots. We also performed a simplified implementaas performance basethe precision of the two first expressing their capacity to give the same senses/Wikipedia articles initially attached to the text fragments anchors. Throughout our experiments, we le annotating every availaoptimal results up depicting the e between two Wikipedia entities for similar optimization technique was also improved comparatively with TAGME, converging faster near %, but the time complexity of this method was by far larger than both our first method and TAGME. The rel and %, optimal e improvement in both accuracy and training set size for fast convergence to optimal results was remarkable. This fact was Method 2 creates a high article of an annotation in the top of ) selection, . In conclusion, we note the necessity of weightformula has proven.

Finally we evaluated our iterative method's precision, (figure 2), varying the maximum iterations parameter. The optimal results for this method's accuracy exceeded 90.71%. We have not developed an optimization technique for the current method due to its numerous parameters, thus its complexity. We applied common precision and recall metrics for the overall evaluation including the pruning step (figure 3), yielding exceptional conclusions for the quality of pruning.

To conclude, method 3 yields less accurate results from the previous methods optimal results, when annotating all available anchors in a text fragment. Yet the optimal accuracy of the current method, is to be explored in future but consist a complex task. In very low pruning recall conditions, accuracy exceeds 98%, with ambiguous anchors contributing as well in these high results. Method 3 accuracy outperforms limited training versions of methods 1, but still remains less accurate than method 2. Method 3 balances speed with accuracy, with the ability to adjust this equilibrium, yet with no expectancy of outperforming method 1 in speed, or method 2 in accuracy.

Figure 3 5 Conclusion

In this paper we propose three novel methods. The first method, revealed the potentiality for further improvement of the relatedness measures between Wikipedia entities. Our second method yields very accurate results, yet requires substantial resources and time. The τ-optimization algorithm we developed for training our two first methods, allowed exploring the optimal upper bound of their performance, and the behavior after training with various training set sizes. The cornerstone of this contribution is a novel iterative approach of the Wikification task, not explored by previous work. This technique, achieves convergence through a series of iterative steps, each of which aims at improving approximation of the optimal annotations.

The development of a better relatedness measure for Wikipedia entities, based on both incoming links, and contextual content is included among our future plans. We also aim at further improving, fine tuning our methods performance and engineering their large scale deployment as an online service, after reserving the necessary resources.

Acknowledgments

This research has been co-financed by the European Union (European Social Fund -ESF) and Greek national funds through the Operational Program "Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF) -Research Funding Program: Thales. Investing in knowledge society through the European Social Fund"


Konstantinos Moustakas
email: moustakas@upatras.gr

Six Degrees of Freedom Implicit Haptic Rendering

Keywords: Haptic rendering, implicit representation, six degrees of freedom

This paper introduces a six degrees of freedom haptic rendering scheme based on an implicit support plane mapping representation of the object geometries. The proposed scheme enables, under specific assumptions, the analytical reconstruction of the rigid 3D object's surface, using the equations of the support planes and their respective distance map. As a direct consequence, the problem of calculating the force feedback can be analytically solved using only information about the 3D object's spatial transformation and position of the haptic probe. Several haptic effects are derived by the proposed mesh-free haptic rendering formulation. Experimental evaluation and computational complexity analysis demonstrates that the proposed approach can reduce significantly the computational cost when compared to existing methods.

Introduction

Human perception combines information of various sensors, including visual, aural, haptic, olfactory, in order to perceive the environment. Virtual reality applications aim to immerse the user into a virtual environment by providing artificial input to its interaction sensors (i.e., eyes, ears, hands). The visual and aural inputs are the most important factors in human-computer interaction (HCI). However, virtual reality applications will remain far from being realistic without providing to the user the sense of touch. The use of haptics augments the standard audiovisual HCI by offering to the user an alternative way of interaction with the virtual environment [START_REF] Burdea | Virtual Reality Technology[END_REF]. However, haptic interaction involves complex and computationally intensive processes, like collision detection or distance calculation [START_REF] Moustakas | Haptic rendering of visual data for the visually impaired[END_REF], that place significant barriers in the generation of accurate and high fidelity force feedback.

Seen from a computational perspective, haptic rendering can be decomposed in two different but heavily interrelated processes, namely collision detection and force calculation. Initially, collisions have to be identified and localized and then the resulting force feedback has to be estimated so as to accurately render the force that will be fed back to the user using specific assumptions on the physical model involved. Concerning haptic rendering research can be divided into three main categories [START_REF] Lin | Haptic rendering: Foundations, algorithms and applications[END_REF]: Machine Haptics, Human Haptics and Computer Haptics [START_REF] Srinivasan | Haptics in virtual environments: Taxonomy, research status and challenges[END_REF]. Machine Haptics is related to the design of haptic devices and interfaces, while Human Haptics is devoted to the study of the human perceptual abilities related to the sense of touch. Computer Haptics, or alternatively haptic rendering, studies the artificial generation and rendering of haptic stimuli for the human user [START_REF] Moustakas | Sq-map: Efficient layered collision detection and haptic rendering[END_REF], [START_REF] Moustakas | Haptic rendering of visual data for the visually impaired[END_REF]. It should be mentioned that the proposed framework takes into account recent research on human haptics, while it provides mathematical tools targeting mainly the area of computer haptics.

The simplest haptic rendering approaches focus on the interaction with the virtual environment using a single point. Many approaches have been proposed so far both for polygonal, non-polygonal models, or even for the artificial generation of surface effects like stiffness, texture or friction, [START_REF] Laycock | A survey of haptic rendering techniques[END_REF]. The assumption, however, of a single interaction point limits the realism of haptic interaction since it is contradictory to the rendering of more complex effects like torque. On contrary, multipoint, or object based haptic rendering approaches use a particular virtual object to interact with the environment and therefore, besides the position of the object, its orientation becomes critical for the rendering of torques. Apart from techniques for polygonal and non-polygonal models [START_REF] Laycock | A survey of haptic rendering techniques[END_REF], voxel based approaches for haptic rendering including volumetric haptic rendering schemes [START_REF] Palmerius | Haptic rendering of dynamic volumetric data[END_REF] have lately emerged. Additionally, research has also tackled with partial success the problem of haptic rendering of dynamic systems like deformable models and fluids [START_REF] Barbic | Six-dof haptic rendering of contact between geometrically complex reduced deformable models: Haptic demo[END_REF].

In general, with the exception of some approaches related to haptic rendering of distance or force fields, one of the biggest bottlenecks of current schemes is that haptic rendering depends on the fast and accurate resolution of collision queries. The proposed approach aims to widen this bottleneck by providing a free-form implicit haptic rendering scheme based on support plane mappings able to provide a six degrees of freedom haptic feedback. In particular, a 3D object is initially modelled using the associated support plane mappings [START_REF] Van Den Bergen | Collision Detection in Interactive 3D Environments[END_REF]. Then the distance of the object's surface from the support plane is mapped at discrete samples on the plane and stored at a preprocessing step following the same procedure presented in [START_REF] Moustakas | Free-form implicit haptic rendering[END_REF]. During run-time and after collision queries are resolved, estimation of the force feedback can be analytically estimated, while several haptic effects, like friction and texture can be easily derived. This results in constant or linear time complexity haptic rendering based only on the 3D transformation of the associated object and the position of the haptic proxy.

Support plane mappings

Support planes are a well studied subject of computational geometry and have been employed in algorithms for the separation of convex objects [START_REF] Van Den Bergen | Collision Detection in Interactive 3D Environments[END_REF][START_REF] Dobkin | A linear algorithm for determining the separation of convex polyhedra[END_REF][START_REF] Chung | Quick Collision Detection of Polytopes in Virtual Environments[END_REF]. From a geometrical perspective, a support plane E of a 3D convex object O is a plane such that O lies entirely on its negative halfspace H - E . Support planes have become useful in previous algorithms based on the concept of support mappings. A support mapping is a function that maps a vector v to the vertex of vert(O) of object O that is "most" parallel to v [START_REF] Van Den Bergen | Collision Detection in Interactive 3D Environments[END_REF][START_REF] Ericson | Real-Time Collision Detection[END_REF]. As a direct consequence, a support plane can be defined as the plane that passes through s O (v), the support mapping of v, and is parallel to v.

The importance of support planes is intuitively apparent: they provide an explicit way of deciding whether another object could possibly intersect with the one that the support planes refers to. Based on this simple but important feature of support planes, a slightly more generalized formulation has been derived [START_REF] Vogiannou | Enhancing bounding volumes using support plane mappings for collision detection[END_REF] introducing the concept of support plane mappings for collision detection. The approach described in [START_REF] Vogiannou | Enhancing bounding volumes using support plane mappings for collision detection[END_REF] is used in the proposed framework for collision detection.

After collision is detected, the force feedback provided to the user through the haptic device has to be calculated. In the present framework, force feedback is obtained directly from the model adopted for collision detection, thus handling collision detection and haptic rendering in an integrated way, as described in the sequel.

Let the parametric form of the support plane equation S SP (η, ω) be:

S SP (η, ω) =   x 0 + ηu 1 + ωv 1 y 0 + ηu 2 + ωv 2 z 0 + ηu 3 + ωv 3   , ∀η, ω ∈ (1)
where u and v constitute an orthonormal basis of the support plane and (x 0 , y 0 , z 0 ) its origin.

Assuming now a dense discretization of the η, ω space, we can define a discrete distance map of the support plane SP and the underlying manifold mesh surface S mesh , by calculating the distance of each point of SP from S mesh :

D SP (η, ω) = ICD (S SP , S mesh ) ( 2
)
where ICD calculates the distance of every point sample (η, ω) of the support plane SP , alongside the normal direction at point (η, ω), from the mesh S mesh and assigns the corresponding values to the distance map D SP (η, ω). The distance map is used in the sequel to analytically estimate the force feedback.

It should be mentioned that the above procedure results in scalar distance maps that accurately encode the surface if and only if there is a one to one mapping of all surface parts with at least one support plane. If such a mapping does not exist, then vectorial distance maps can be used that include information about the distance of all sections of the ray cast in the normal direction of the support plane to the object mesh.

Point-based (3DoF) haptic rendering

Referring to Figure 1, let point H p be the position of the haptic probe and S mesh represent the local surface of the object.

Let also S SP represent the distance of point H p from the support plane, which corresponds to point P M on the SP.

If collision is detected, the absolute value of the force fed onto the haptic device is obtained using a spring model as illustrated in Figure 1. In particular:

F = k • |S SP -D SP (P M )| ( 3
)
where k is the spring constant. D SP (P M ) is the distance of point P M from the mesh and is stored in the distance map of the support plane. Notice that the term |S SP -D SP (P M )| is an approximation of the actual distance of H p from the mesh that becomes more accurate if the support plane surface approximates well the mesh.

The direction of the force should in general be perpendicular to the local area, where collision is detected. An obvious solution to the evaluation of the direction of this force would be to detect the surface element (i.e. triangle) where the collision occurred and to provide the feedback perpendicularly to it. This approach is not only computationally intensive, but also results in non-realistic noncontinuous forces at the surface element boundaries. In the present framework the analytical approximation of the mesh surface is used utilizing the already obtained SP approximation and the distance map. Based on this approximation the normal to the object's surface can be approximated rapidly with high accuracy. In particular, assuming that (η, ω) form an orthonormal basis of the plane then if D SP (η, ω) is the scalar function of the distance map on the support plane, as previously described, the surface S mesh of the modelled object can be approximated by equation (4) (Figure 1):

S mesh (η, ω) = S SP (η, ω) -D SP (η, ω) n SP (4)
where S SP is the surface of the support plane, D SP the associated distance map and n SP its normal vector that can be easily evaluated through n SP = u × v. Now the calculation of the force feedback demands the evaluation of the normal vector n S on the object's surface. That is obtained through equation [START_REF] Moustakas | Haptic rendering of visual data for the visually impaired[END_REF]. In the following the brackets (η, ω) will be omitted for the sake of simplicity.

n S = ∂S mesh ∂η × ∂S mesh ∂ω (5)
where

∂S mesh ∂η = ∂S SP ∂η - ∂D SP ∂η n SP -D SP ∂n SP ∂η (6)
Since n SP is constant over SP, equation ( 6) becomes:

∂S mesh ∂η = u - ∂D SP ∂η n SP (7)
A similar formula can be extracted for ∂S mesh ∂ω :

∂S mesh ∂ω = v - ∂D SP ∂ω n SP (8)
All above terms can be computed analytically, except from ∂D SP ∂η and ∂D SP ∂ω that are computed numerically.

Substituting now equations ( 4), ( 6), ( 7), [START_REF] Barbic | Six-dof haptic rendering of contact between geometrically complex reduced deformable models: Haptic demo[END_REF] in equation ( 5) the normal direction n S can be obtained.

Since, the direction of the normal along the surface of the modelled object is obtained using equation ( 5), the resulting force feedback is calculated through:

F h = k |S SP -D SP (P M )| n S n S (9)

6DoF haptic rendering

Let us now assume that the haptic interaction point is actually a haptic interaction object that is also modelled using support plane mappings and distance maps.

Referring to Figure 2 that for simplicity depicts the 2D case, let M 1 and M 2 be the mesh areas of two different objects possibly involved in collision and S 1 and S 2 their respective support planes. Based on the approach presented in [START_REF] Vogiannou | Enhancing bounding volumes using support plane mappings for collision detection[END_REF] it can be decided whether M 1 and M 2 are probably involved in collisions using their SPMs S 1 and S 2 . Now the question is: Is it possible to identify the collisions and calculate 6DoF force feedback without entering the computationally intensive narrow-phase [START_REF] Vogiannou | Enhancing bounding volumes using support plane mappings for collision detection[END_REF] of the collision detection algorithm?

This question reduces to the problem of calculating the impact volume (3D case) or impact surface (2D case) S as depicted in Figure 2. Let us consider for sake of simplicity and without loss of generality the 2D case. If d is a material density function then the mass corresponding to the collision area S can be described by the following formula.

V = S ddS (10)
Now referring again to Figure 2 let W be a sampling point of the support plane S 2 and W 2 the point of mesh M 2 that can be trivially reconstructed by projecting W along side the normal direction n 2 as far as the distance map Fig. 2. Distance calculation using distance maps over support planes D SP (W ) dictates. Now, W 2 can be also trivially projected on the support plane S 1 and retrieve the corresponding distance value D 1 (M 1 ) of S 1 from the mesh M 1 . Now, omitting the notation (η, ω) in the following for simplicity and assuming the function f is defined as follows:

f = D 1 (W 2 ) -D 1 (M 1 ) (11)
where D 1 (W 2 ) is the distance of point W 2 from S 1 and D 1 (M 1 ) the corresponding distance of S 1 from mesh M 1 . Now, collision can be reported only for cases where function f is positive. Let now C 1 and C 2 (Figure 2) denote the projection of the colliding volume on the support planes S 1 and S 2 respectively. It is obvious that C 1 is the support set of positive values of function f . Then, equation ( 10) can be transformed as follows:

V = t1∈C1 d • f dt 1 = t2∈C2 d • f • cos (n 1 n 2 ) dt 2 (12)
A similar analysis can be followed for the three dimensions and the corresponding formula is:

V = η,ω∈C2 d • f • cos (n 1 n 2 ) dηdω (13)
Now since the volume of the penetrating object part and its centre of inertia is known, 6DoF haptic feedback can be easily estimated. In the context of the proposed framework a single force is estimated for the penetrating part of the haptic interaction object that is applied on the gravity centre of the colliding volume. Other approaches for approximate 6DoF haptic rendering like the one presented in [START_REF] Mcneely | Six degree-of-freedom haptic rendering using voxel sampling[END_REF] can be also implemented.

Haptic effects

The analytical estimation of the force feedback based only on the object 3D transformation, the probe position and the distance maps, provides the opportunity to develop closed form solutions for the rendering of physics-based or symbolic force effects; the following sections indicatively describe some of them.

By applying a local smoothing operation to the distance map, the resulting force feedback is smooth in the areas around the edges, without being overrounded as is the case with the force shading method [START_REF] Ruspini | The haptic display of complex graphical environments[END_REF]. A typical example of distance map preprocessing so as to achieve force smoothing using a Gaussian kernel is given by the following equation:

D SP (η, ω) = D SP (η, ω) * G σ (η, ω) (14)
where G σ is a 2D Gaussian kernel and " * " denotes convolution. It is evident that different smoothing operators can be easily applied. A very useful operator that can be implemented is the force smoothing only in areas that are not smooth due to the finite tessellation (sampling) and not in object and surface boundaries, following the popular in computer graphics "crease angle" concept. A haptic "crease angle" rendering can be easily performed by applying anisotropic diffusion or even an edge-preserving smoothing operator on the distance map. Similarly using the proposed framework for haptic rendering, haptic texture can be also simulated easily by applying appropriate transformations on the acquired distance map. An example for simulating surface roughness is provided below, where Gaussian noise is added on the distance map. No computational cost is added, since the procedures for calculating the force direction are not altered due to the existence of haptic texture. The only difference lies in the evaluation of the magnitude of F texture , which now yields from:

F texture = k |S SP -(D SP (P M ) + n g )| n S n S (15)
where n g denotes the gaussian noise.

Complexity and experimental results

In the following an analysis of the computational complexity of the proposed scheme in comparison to the typical state-of-the-art mesh-based haptic rendering scheme is discussed. Moreover, even if an experimental analysis of the proposed support plane mapping based haptic rendering approach, in terms of timings for simulation benchmarks would not be fair for the state-of-the-art approaches, since it would encode the superiority of SPM based collision detection and would not directly highlight the proposed haptic rendering approach, two experiments are presented where the proposed haptic rendering scheme is compared to the state-of-the-art mesh-based haptic rendering in terms of timings in computationally intensive surface sliding experiments.

After collision is reported, a typical force feedback calculation scheme would need to identify the colliding triangle of the involved 3D object in O(n) time, where n is the number of triangles, or in O(logn) time if bounding volume hierarchies are used. Then the force can be calculated in constant O(1) time. In order to avoid force discontinuities, for example force shading, and if there is no adjacency information then the local neighborhood of the colliding triangle can be found again in O(n) time, where n is the number of triangles, or in O(logn) time if bounding volume hierarchies are used. Finally, the mesh-based haptic rendering scheme has no additional memory requirements per se. On the other hand, concerning the proposed free-form implicit haptic rendering scheme, after a collision is detected, the resulting force feedback can be calculated in constant time O(1) using equation [START_REF] Van Den Bergen | Collision Detection in Interactive 3D Environments[END_REF]. In order to avoid depth discontinuities the distance map can be smoothed, in an image processing sense, at a preprocessing phase. Even if this step is performed during run-time it would take O(k) time, where k is the local smoothing region or the filtering kernel window. Concerning the 6DoF case, it can be considered requiring to execute m times the 3DoF procedure, where m is the number of support plane sampling points processed. On the other hand the proposed scheme has O(m • s) memory requirements, where m is the number of support planes and s the number of samples per support plane. Taking now into account that the more support planes are used the smaller their size and the less samples are necessary for a specific sampling density we can safely assume that the memory requirements are linear to the total number of samples that depends on the sampling density used.

Table 1 summarizes the computational complexity analysis of the proposed free-form haptic rendering scheme, when compared to the mesh-based approach. Concerning the quantitative results, interaction with two objects was considered, namely the Galeon and Teeth models of 4698 and 23000 triangles respectively. The objects are illustrated in Figure 3 .

Moreover, CHAI-3D was used for interfacing with the haptic devices [START_REF] Conti | Chai 3d: An open-source library for the rapid development of haptic scenes[END_REF]. Moreover, the force estimation algorithms were applied on a predefined trajectory of the haptic probe, so as to assure fair comparison. In particular, initially the trajectory of the haptic probe in the 3D space has been recorded while being in sliding motion over the objects' surface. Then this trajectory has been used as input for both algorithms so as to extract the timings mentioned below. Table 2 presents the mean timings and their standard deviation of the force estima- It should be emphasized that the above timings need to be taken into account under the exact experimental setting. In particular, concerning the proposed approach 1000 support planes were used for the case of the Galeon and 1500 for the case of the Teeth model. Distances are estimated for all support planes and forces are calculated for the closer one. This procedure, could be optimized by partitioning the space in a preprocessing step and knowing beforehand to which support plane, each point in space "belongs to", thus reducing the search from O(n) to O(logn). Moreover, concerning the mesh-based approach force shading has been also implemented.

Conclusions

The proposed framework extends the support plane mapping concept and the corresponding collision detection scheme to direct free-form implicit haptic rendering. An analytical scheme to calculate both 3DoF and 6DoF force feedback is proposed. Its most significant property is that there is no need to enter the narrow-phase of the collision detection pipeline so as to identify the colliding triangles and subsequently estimate the reaction forces, on contrary it requires only the distance maps stored on the support planes. It is evident that the proposed scheme reduces significantly the computational cost in the performed simulations. This significant gain comes at an expense of two limitations. Firstly, special care has to be taken at the preprocessing step so that the models are well approximated using the support planes and the distance maps. For example if the objects demonstrate large complex concavities the use of vectorial distance maps is inevitable. Secondly, the proposed scheme cannot be, in its current form,

Fig. 1 .

 1 Fig. 1. Distance calculation using distance maps over support planes

Fig. 3 .

 3 Fig. 3. Galeon model, 4698 triangles. Teeth model, 23000 triangles

Table 1 .

 1 Computational complexity comparison

	Process	Mesh-based	Free-form
	Force	O(n) or O(logn)	O(1)
	6DoF force O(mn) or O(mlogn) O(m)
	Smoothing O(n) or O(logn)	O(1)
	Memory	-	O(m • s)

Table 2 .

 2 Galeon and Teeth models: Interaction timings

		Galeon		Teeth
	Process Mean time (ms) σ Mean time (ms) σ
	Mesh-based	1.2	0.22	4.2	0.82
	Free-form	0.028	0.006	0.036	0.005

directly applied to deformable models. An extension to piecewise or free-form deformable models, where deformations can be analytically expressed seems possible and remains a direction for future work.