Valé Rie Renaudin
email: valerie.renaudin@ifsttar.fr

Melania Susi
email: msusi@ucalgary.ca

Gé Rard Lachapelle
email: gerard.lachapelle@ucalgary.ca

Step Length Estimation Using Handheld Inertial Sensors

Keywords: biomechanics, dead reckoning, pedestrian navigation, step length, handheld devices, IMU

In this paper a novel step length model using a handheld Micro Electrical Mechanical System (MEMS) is presented. It combines the user's step frequency and height with a set of three parameters for estimating step length. The model has been developed and trained using 12 different subjects: six men and six women. For reliable estimation of the step frequency with a handheld device, the frequency content of the handheld sensor's signal is extracted by applying the Short Time Fourier Transform (STFT) independently from the step detection process. The relationship between step and hand frequencies is analyzed for different hand's motions and sensor carrying modes. For this purpose, the frequency content of synchronized signals collected with two sensors placed in the hand and on the foot of a pedestrian has been extracted. Performance of the proposed step length model is assessed with several field tests involving 10 test subjects different from the above 12. The percentages of error over the travelled distance using universal parameters and a set of parameters calibrated for each subject are compared. The fitted solutions show an error between 2.5 and 5% of the travelled distance, which is comparable with that achieved by models proposed in the literature for body fixed sensors only.

Introduction

Recent progress in Micro Electro Mechanical Systems (MEMS) technology is stimulating their use in different domains including pedestrian navigation, location based services (LBS), safety and healthcare services. Because they are already embedded in many electronic devices, and particularly in smart phones, it is now possible to use these low cost sensors for developing new services, pedestrian navigation being one of the most attractive for the consumer grade market. Indeed MEMS accelerometers and gyroscopes provide autonomous solutions for tracking pedestrians in different types of environments, thereby overcoming the limitations of Global Navigation Satellite Systems (GNSS) in challenging contexts, such as indoors or urban canyons, where satellite signals are blocked or strongly attenuated by man-made infrastructures. Beyond using pedestrian navigation for commercial applications, the ability to continuously track people anywhere can answer safety issues such as monitoring vulnerable patients. This is, for example, the case for subject affected by a cognitive function deficit. Helping elderly people suffering from dementia or Alzheimer's is a good illustration. Indeed the main effects of these pathologies are memory loss and attention deficit inducing difficult orientation and wandering. These kinds of patients require continuous tracking and monitoring to provide prompt assistance in case of necessity or to alert third parties when they wander beyond a specific radius considered as a -safe‖ zone. However, MEMS sensors cannot provide an accurate self-contained system mainly due to the errors inherent to their low-cost nature, namely drift and biases affecting their signals [START_REF] Titterton | Strapdown Inertial Navigation Technology[END_REF]. Frequent external sources of information, for example GNSS updates, are needed for mitigating the sensors errors.

When GNSS aiding is not available alternative approaches must be adopted. For pedestrian navigation, the characteristics of human gait can be exploited in Pedestrian Dead Reckoning (PDR) algorithms instead of double integrating the inertial data, which is implemented in the strap-down method. Indeed the latter approach is not suitable for low-cost sensor based applications, since the double integration increases the noise component proportionally to the operational time even if the pedestrian is not moving. Instead in a PDR approach, the estimation of the current pedestrian's position results from the displacement of the user, i.e., linear walked distance and walking direction, since the last known position estimate. This recursive process is related to the effective motion of the user.

The computation of the user's linear displacement generally consists of two parts: first detecting the user's steps and second evaluating their length. Estimating a pedestrian's step length is a challenging task that can be performed following different approaches, which strongly depend on the sensor's location. The majority of existing algorithms assumes that the sensor is rigidly attached to the user's body either on the foot, close to the Centre Of Mass (COM), e.g., along the backbone, or distributed on the leg [START_REF] Renaudin | Indoor navigation of emergency agents European[END_REF][START_REF] Tien | Results of using a wireless inertial measuring system to quantify gait motions in control subjects[END_REF][START_REF] Miyazaki | Long-term unrestrained measurement of stride length and walking velocity utilizing a piezoelectric gyroscope[END_REF][START_REF] Jahn | Comparison and Evaluation of Acceleration Based Step Length Estimators for Handheld Devices[END_REF][START_REF] Alvarez | Comparison of Step Length Estimators from Wearable Accelerometer Devices[END_REF]. These locations are particularly suitable for navigation purposes since the inertial force experienced by the sensor is directly linked to the gait cycle. Using body fixed sensors, two main categories of step length models can be identified in the literature: biomechanical and parametric models. In general biomechanical models assume that the sensor is located on the user's COM and model the user's leg as an inverted pendulum [START_REF] Jahn | Comparison and Evaluation of Acceleration Based Step Length Estimators for Handheld Devices[END_REF][START_REF] Alvarez | Comparison of Step Length Estimators from Wearable Accelerometer Devices[END_REF]. A simple geometric relationship between the COM's vertical displacement and the step length is then applied. Models based on other geometric considerations are also proposed in [START_REF] Zijlstra | Displacement of the pelvis during human walking: experimental data and model predictions[END_REF][START_REF] Kim | A step, stride and heading determination for the pedestrian navigation system[END_REF]. Parametric models use the step frequency and the accelerometers variance, either combined or independently, to estimate the step length [START_REF] Shin | Adaptive Step Length Estimation Algorithm Using Low-Cost MEMS Inertial Sensors[END_REF][START_REF] Sun | Pedestrian dead reckoning based on activity recognition and stride assessment[END_REF]. Again the sensor is either mounted on the belt or on the foot but body fixed locations are not suitable for many applications. As explained, MEMS are often already embedded in unobtrusive portable devices, e.g., smart phones or personal digital assistants, which are usually carried in hands or kept in bags and therefore are -non-body fixed‖. Published work on using non-body fixed sensors for pedestrian navigation is however often constraining the sensor's location to emplacements where the device is relatively stable while the user is walking. For example the device is carried in the user's trouser pocket [START_REF] Steinhoff | Dead Reckoning from the Pocket-An Experimental Study[END_REF] or constrained to specific locations like close to the ear while phoning or pointing toward the walking direction [START_REF] Skog | Zero-velocity detection-An algorithm evaluation[END_REF]. The reason is that in these scenarios, the IMU (Inertial Measurement Unit) signal patterns of the device are closer to the ones produced by body fixed sensors and subsequently similar approaches can be adopted.

When the sensor is handheld without any constraint, the situation becomes much more complex adding many new issues that require specific processing. For example, since the hand undergoes many motions which don't reflect the user's displacement, they have to be identified and classified as parasite in order to avoid wrong propagation of the user's position. In hand, not only is the sensor's orientation unknown a priori but also it can vary very quickly due to fast hand motions. In fact very few studies target the handheld sensor case and in general only the case of sensors held in the user's phoning or texting hand is considered [START_REF] Gusenbauer | Self-Contained Indoor Positioning on Off-The-Shelf Mobile Devices[END_REF]. Indeed in this context, the sensor is mainly experiencing the inertial force produced by the global motion of the user, which is similar to the body fixed case. Conversely, the cases of the sensor held in the swinging hand and when the sensor's placement changes while the user is moving are omitted. In [START_REF] Ayub | Sensor placement modes for smartphone based pedestrian dead reckoning[END_REF], different sensor carrying modes are examined, including carrying the sensor in the swinging hand, but only traditional techniques, designed for body fixed sensors, are adopted. When the above techniques are applied to handheld smart phones, they produce lower performance than the ones obtained with body fixed sensors.

Facing the identified limitations of existing techniques in the context of autonomous indoor navigation based on smart phones, a dedicated and extensive analysis of the hand case has been performed herein. Its results are presented in this paper and lead to the development of a handheld based step length model. Algorithms are proposed for estimating the step length of pedestrians walking on a flat ground using handheld devices without constraining the sensor's carrying mode. The proposed step length model combines the user's step frequency and height. Step frequency evaluation is performed directly in the frequency domain and independently from the step detection process. In order to adapt the model to the handheld case, the relationship between step frequency and hand frequency is deeply investigated. Performance of the proposed model is assessed in the position domain by combining the step length model with a step detection algorithm presented in [START_REF] Susi | Motion mode recognition and step detection algorithms for mobile phone user[END_REF]. The assessment part, performed with 10 test subjects, shows that the handheld step length model achieves comparable performances as the ones obtained in the literature but with body fixed sensors only.

The structure of the paper is the following. In Section 2, the signal model is introduced and the signal preprocessing phase is illustrated. In Section 3, the analysis of human gait using handheld devices is described. Then, in Section 4, the proposed step length model is presented with a description of a novel technique used to extract the user's step frequency from the user's hand frequency. Section 5 deals with the assessment of the proposed algorithm with 10 test subjects. Finally, Section 6 draws conclusions.

Signal Model and Pre-Processing

In this paper step length estimation is performed using a six-degree of freedom (6DoF) IMU. It comprises a tri-axis gyroscope and accelerometer that sense angular rates and accelerations of the body frame. The sensor's digital output is modeled as a six-dimensional vector given by the sum of the response to the sensed inertial force and a noise term [START_REF] Skog | Zero-velocity detection-An algorithm evaluation[END_REF]:

              aa n n n n n nn                          sa s s    (1)
where:

 n  ℕ,
  n   .
Since the frequency content of the accelerations and angular rates induced by human gait are below 15 Hz [START_REF] Mathie | Monitoring and Interpreting Human Movement Patterns Using a Triaxial Accelerometer[END_REF], the components in Equation ( 1) are low pass filtered using a zero-lag 10th order Butterworth filter with a 15 Hz cut-off frequency. It is worth underlining that although fast hand motions produce frequencies larger than 15 Hz, these events are not of interest for pedestrian navigation because they are not related to the user's global locomotion and consequently can be removed. The filtered components are indicated with   When the sensor is not body fixed, the orientation of the handheld device is a priori unknown, which implies that the signal processing has to be immune to any change of the sensor's orientation. This is even more important with a handheld device as the rapid motions characteristic of the hand render the task of continuously estimating the device's orientation very challenging. Therefore, instead of working on the individual vector elements, the norm of the filtered components is considered:

      a rms rms rms sn n sn       s (2)
where .

Finally, because the zero frequency component (DC), which is present in the signals, can mask information related to the different motion modes experienced by the user and reduce the effectiveness of the signal frequency analysis, it has been removed from the signal of interest: [START_REF] Tien | Results of using a wireless inertial measuring system to quantify gait motions in control subjects[END_REF] where the second term is the signal's mean evaluated using a moving average filter. The selection of the window's length (L) is described later.

  n  0 1 0 1 ] [ ] [ ] [ rm L rms s rms l s n l s n s L n      

Gait Analysis for Step Length Estimation Using Handheld Devices

Before evaluating the step length, step events must be identified. When the sensor is foot mounted, step events can be detected using periods of zero velocity or zero angular rate corresponding to the stance phase, i.e., the period during which the foot is flat on the ground [START_REF] Skog | Zero-velocity detection-An algorithm evaluation[END_REF]. With handheld devices, these periods are not present. In addition, when the sensor is not body fixed, signal patterns can suddenly vary following a change of the carrying mode or the user's motion. Therefore knowing the nature of the user's motion and the device's carrying mode is critical for designing adaptive algorithms and improving the robustness. To cope with the complexity of the hand case, the process for estimating the user's travelled distance has been divided in three phases: motion mode recognition, step detection and step length estimation. Motion mode recognition and step event detection are performed applying the approach described in details in [START_REF] Susi | Motion mode recognition and step detection algorithms for mobile phone user[END_REF]. This approach is now summarized in Sections 3.1 to 3.5, in order to provide a complete description of the novel step length estimation procedure based on a handheld device.

Motion Mode Recognition

The recognition of the motions of the subject and his/her hand is considered as a classification problem. Six motion modes that are typical for portable device users have been identified, namely:

(1) Static: the subject is not moving even if a slight motion occurs.

(2) Walking with a swinging hand: the user is walking holding the portable device in the hand.

(3) Walking with a texting hand: the user is walking while typing or reading a message/instruction on the screen's device. (4) Walking with phoning hand: the user is using the device to make or receive a phone call.

(5) Walking with the mobile device in the handbag. ( 6) Irregular motion: this class includes all motions that are not related to a real displacement change of the user.

It has been found that texting, phoning and carrying the device can be considered as a unique class for the purpose to optimize the step detection process [START_REF] Susi | Motion mode recognition and step detection algorithms for mobile phone user[END_REF]. During these three motions, the inertial force sensed by the IMU is primarily due to the lower part of the user's body and the sensor signals present patterns that are close to those recorded with body fixed sensors. Consequently, the number of classes describing the motion modes of interest is reduced to four. The classification process is performed by extracting a set of features from the raw sensor data. These features are computed by dividing the data in windows of 256 samples, corresponding to 2.56 seconds, with a 50% overlap. The window size, indicated as L in (3), is selected to be small enough to catch any fast motion's change and large enough to include one complete gait cycle. In addition, a window size 2 n N  allows the computation of the Fast Fourier Transform (FFT) used for the analysis in the frequency domain. It is worth mentioning that the proposed FFT analysis is applied on a window that is small enough to target real time applications since the induced lag is around half of the window size, namely 1.28 seconds. Real time functionalities have not yet been implemented. The following features, both in time and frequency domains, are extracted to identify the motion mode: the energies, the variances and dominant frequencies of the IMU signal.

Time Domain Features: Energy Related Features and Signal's Variance

The energy features enable distinguishing between activities of low and high intensity. The energy is evaluated by computing the norm of the accelerometer and the gyroscope's pre-processed measurements, adding and normalizing them over the analysis window as:

  0 2 1 0 1 a rms N a rms s n E s n N     (4)   0 2 1 0 1 rms N rms s n E s n N       (5)
where N is the length of the analysis window The swinging mode experiences much higher amplitudes of angular rate and acceleration energies than the other states, which are the texting, phoning and bag carrying motion modes. Consequently, the variance of the gyroscope signal is also used to distinguish the above states. To refine the classes' characterization, the variances of the gyroscope and accelerometer signals are also evaluated. This statistical measurement is defined as the average of the squared differences from the mean: This allows the swinging mode to be identified as the variances of both gyroscope and accelerometer signals are larger than in the other cases. This feature is also used to recognize irregular motions characterized by a sudden increase of the variance for both inertial data but without observing any signal periodicity.

   

Frequency Domain Features: Dominant Frequencies

The study of human gait has shown that walking is characterized by a fundamental pattern that is not subjected to inter and intra individual variations and that is produced by the cyclic repetition of the stance and the swinging phases of the foot [START_REF] Rose | Human Walking[END_REF][START_REF] Sabatini | Assessment of walking features from foot inertial sensing[END_REF]. By analyzing the accelerations in the frequency domain, the human gait periodicity can be captured. Since the signal produced by the walk is not stationary, i.e., the signal's statistical properties change over time, the frequency analysis is performed using the Short Time Fourier Transform (STFT) [START_REF] Cohen | Time Frequency Analysis: Theory and Applications[END_REF]. This signal processing technique is based on the fragmentation of the input signal in short temporal windows where the signal is assumed to be stationary. The STFT approach, despite its poor time-frequency localization properties, has been selected for its low computational cost and for its non-parametric nature. The spectrogram, a time varying spectral representation of the signal, is then obtained by squaring the absolute value of the STFT. In Figure 1, the spectrogram of the accelerometer signal is reported for a walking subject with the sensor on the foot. Three dominant frequencies, i.e., maximal frequencies in the spectrogram produced by the three main temporal periodicities, are clearly visible. As shown in Figure 2, when the sensor is held in hand, three dominant frequencies are also identifiable in the spectrogram of the accelerometer signal for a subject walking alternatively with the sensor in the swinging and in the texting modes. As shown in [START_REF] Susi | Motion mode recognition and step detection algorithms for mobile phone user[END_REF][START_REF] Susi | Quasi-Static Detection from MEMS Handheld Devices[END_REF], the first two dominant frequencies can be used to distinguish irregular motion and walking mode. Similar analysis can be conducted with the handheld gyroscope signal to distinguish between different motion modes. Indeed a user walking with the sensor in the swinging hand induces peaks in the spectrogram of the gyroscope signal, due to the periodic rotation of the arm during its swinging phase. Conversely, when the sensor is in the texting mode, phoning hand mode or in the user's bag, dominant peaks are not produced by the gyroscope signal.

Sensors 2012, 12 8514

Motion Mode Decision Tree

As shown in Figure 3, all above features are integrated in a multivariate decision tree classifier. In Figure 3 the three phases, namely motion mode recognition, step detection and step length estimation necessary to estimate the user's travelled distance, are also represented. The tree classifier selects one state among all possible motion states represented by the tree leaves. The selection is performed at each node through several tests using the features previously described. The decision tree initially differentiates static and dynamic activity using the energies and variances of MEMS signals. If the activity is classified as dynamic, the periodicity of the accelerations is analyzed. The latter reflects the periodicities of human gait making it possible to distinguish normal walking and irregular motions. Irregular motion is generally characterized by very high values of the variance in a short period of time, e.g., when the user is pulling out a device from a bag without moving. Finally, the features described in Sections 3.2-3.3 enable distinguishing the states of swinging and texting, phoning and carrying the bag. The classifier achieves high performance independently from the considered subject because the selected features are characterized by a high inter-class distance. More details about the classification process can be found in [START_REF] Susi | Motion mode recognition and step detection algorithms for mobile phone user[END_REF].

Step Identification

Once the motion state is determined, an adaptive step detection algorithm is applied. If the IMU is held in a swinging hand, the experienced inertial force results from the gross motion of the subject plus the swinging motion. However biomechanical studies have shown that a synchronization exists between the swinging of the arm and the foot motion during walking. Specifically, it has been observed that during normal walking, the arm swing is produced for decreasing the reaction momentum about the vertical axis of the foot [START_REF] Park | Synthesis of natural arm swing motion in human bipedal walking[END_REF]. When the left foot is in the stance phase, a positive torque, produced by the arm swing, allows the advancement of the right leg. This relation, which is valid in the absence of particular pathologies, allows the identification of step events from the swing of the user's arm. The step information is extracted by using angular rates because the periodic rotation of the arm produces a sinusoidal pattern in the gyroscopes' signals. Subsequently, peaks in the gyroscopes' signals are determined and the up and forward hand's motions are detected along with the synchronized swing phase of the foot [START_REF] Susi | Motion mode recognition and step detection algorithms for mobile phone user[END_REF].

When the sensor is placed in the user's texting/phoning hand or carried in the bag/trouser pocket, it is mainly reflecting the general motion of the user. Even in these cases, periods of zero velocity cannot be determined in the accelerometer signal. However the repetition of swing phases produces negative and positive peaks in the accelerometer signal pattern. These peaks are used to mark step occurrences. Signal peak detection is performed by recognizing a local maximum or minimum within the sliding window and by exploiting algorithms based on adaptive threshold. Thanks to the adaptive feature, the peak detection algorithm is independent from the level of the signal energy and therefore the variations that the hand could undergo. To minimize the probability of false peak detection, a dedicated preprocessing phase of IMU signals has been added. Inertial signals are low-pass filtered using a 10th order Butterworth filter with a 3-Hz cut-off frequency, which produces an undistorted signal for detecting the fundamental frequency produced by the step events. The cut-off frequency value is selected considering the typical range of step frequencies experienced by a pedestrian walking at a normal speed [START_REF] Levi | Dead Reckoning Navigational System Using Accelerometer to Measure Foot Impacts[END_REF]. Finally the algorithm adopts the mean value over a sliding window as the threshold for detecting peaks. If a sample in the window gives a larger value than the adaptive threshold, a peak is identified and a step is detected.

Step Length Model

Once a step is detected, the next stage for tracking pedestrians is to determine their position using the estimated step length. With body fixed sensors, it has been shown in literature that a linear relationship between step length and step frequency exists [START_REF] Shin | Adaptive Step Length Estimation Algorithm Using Low-Cost MEMS Inertial Sensors[END_REF][START_REF] Levi | Dead Reckoning Navigational System Using Accelerometer to Measure Foot Impacts[END_REF]. This can be intuitively understood and expanded to handheld devices. If a pedestrian walks faster, both step's length and step frequency will increase. Biomechanical studies have shown that in general, the user's step length is proportional to the length of the user's leg and subsequently to the user's height [START_REF] Rose | Human Walking[END_REF]. Starting from these results, a new step length model for handheld devices has been designed and empirically validated. It combines the step frequency and user's height.

The best linear relationship between frequency and step length was found to be the one weighted by the user's height and is given by:

+ c (8) ℝ

where h is the user's height, f step is the step frequency and K is a set of parameters. A universal and a calibrated model are proposed. The universal model has been developed for giving a first approximation for any filter that would offer tuning functionalities. As detailed in Section 5 the universal model is based on a set of constants trained using 12 test subjects while the calibrated model tailors the set of constants for each subject.

Recursive Least-Squares (RLS) [START_REF] Teunissen | Adjustment Theory: An Introduction[END_REF] is used to determine the calibrated set of model's parameters. The approach is based on the recursive evaluation of the optimum parameters by minimizing the sum of squared residuals between the true step length and the predicted step length. The universal set of parameters K is used for the initial solution x 0 : [START_REF] Shin | Adaptive Step Length Estimation Algorithm Using Low-Cost MEMS Inertial Sensors[END_REF] with: s comprises the true step lengths for all epochs between 1 and k. H comprises the step frequencies and users' heights for n test subjects over k epochs. The set of fitted parameters is determined when the convergence criteria over is achieved. True step lengths are evaluated following the procedure described in Section 5.1 while predicted step lengths are obtained using the model in Equation (8).

  1 1 1 1 1 0 , , ,1 , , 1 , 1 .

Step Frequency Evaluation

The proposed step length model strongly relies on the extraction of the user's step frequency and its quality. In most of the works, the step frequency is calculated in the time domain by detecting steps and computing the inverse of their duration [START_REF] Kim | A step, stride and heading determination for the pedestrian navigation system[END_REF]. In order to render the algorithms more robust, the step's cadence estimation proposed herein was designed independently from the step detection process. Indeed it is directly evaluated in the frequency domain by computing the FFT of handheld accelerometer signals. This is a critical observation as it ensures that the algorithms constituting the step length estimation process remain uncorrelated.

Since the sensor is located in the hand, the step's frequency cannot be directly extracted from the hand signals. In order to tailor the model to the handheld case, a dedicated analysis was conducted to accurately relate the handheld signal frequencies to the walking cycle. Work concentrated on the coordination between legs, arms and hands.

The examined signals are those of two IMUs, the first one placed in the hand and the second, rigidly fixed on the foot. In Figure 4, the Power Spectral Densities (PSDs) are shown for the accelerometers of both IMUs in the swinging and texting cases. It is observed that the dominant frequency peaks are centered on the same values for both sensors, which shows that the step frequency can be derived even if the sensor is not located on the foot or is not body fixed. Furthermore, the strongest frequency of the accelerometer signal, i.e., the frequency with the maximum power, is not always coupled with the same event of the walking gait cycle. Sometimes it is coupled with the step event (Figure 4) and sometimes, especially for faster speeds, with the stride event. This is further illustrated in Figure 5 where the PSD extracted from the signal of a sensor in the user's swinging hand is reported. The PSD is estimated using the Welch periodogram technique [START_REF] Welch | The use of Fast Fourier Transform for the estimation of the power spectra: A method based on time averaging over short modified periodograms[END_REF]. It is performed by dividing the signal into overlapping blocks and averaging the squared magnitude FFTs (Fast Fourier Transform) of the signal blocks. When compared with the use of the standard periodogram, the main advantage of this method is to minimize the variance, which renders this spectrum estimator unbiased. The spectrum analysis shows that the faster the speed is, the greater the chances are that the strongest frequency and stride event are coupled. Here, the strongest frequency is coupled with stride events. Consequently a binary classifier was designed to identify step frequencies from the extracted hand frequencies. After selecting the strongest frequency, the classifier applies the following decision rules: [START_REF] Steinhoff | Dead Reckoning from the Pocket-An Experimental Study[END_REF] A 1.4 Hz threshold τ has been selected based on the fact that for the normal walking case the range of frequencies for a human step is generally above 1.6 Hz. This assumption has been validated experimentally by analyzing the steps' length of the test subjects. However, further analysis should be performed to investigate how the users' age affects step frequency. In the proposed algorithm, if the detected frequency is the stride frequency, the step frequency can be derived by multiplying the strongest frequency value by two. The step frequencies, extracted after applying the binary classifier, true step lengths and the estimated ones by using the step length model and the universal set of parameters are shown in Figure 6 for a user walking with the sensor located in the swinging hand.

Experimental Assessment

Data Collections

In order to train and assess the performance of the proposed step length model, several field tests were performed. The first type of data collections has been performed to fit the three constants characterizing the set of parameters, indicated as K in Equation ( 8), defining the universal model. The above parameters have been computed using data collected by 12 different subjects: six men and six women between twenty-forty years old. They walked along a 200 m straight line on a parking lot, at three different speeds, for a total of approximately 600 m. True step lengths were measured using a wheel speed sensor whose records were time tagged with GPS time and inertial signals from a foot mounted IMU. As shown in the left side of Figure 7, the wheel speed sensor was controlled by an instructor, whose role was to pace the pedestrian. The test subjects were requested to walk at slow speed (about 0.8 km/h), intermediate speed (about 1.8 km/h) and fast speed (about 4 km/h) with the hand in texting and swinging motion modes. Foot mounted inertial signals were also used to assess the

f f f f f f           
proposed step detection algorithms with a handheld IMU. Detection of the foot stance phases was performed by assessing the acceleration variances. Universal parameters are the outcome of fitting Equation ( 8) to all 12 datasets simultaneously. Figure 8 shows the best fitting of the step length model as a function of the user's height and the product of the user's height with the strongest gait frequency extracted from the handheld device. The true step lengths are depicted in blue. The higher the latter product is the faster the pedestrian walking speed is. For assessing the performance of the proposed model in the position domain, a second experiment was conducted in an open soccer field with different test subjects than the ones who participated to the fitting of the universal step length model. All data collections were performed using a multi-sensor navigation platform, the NavCube [START_REF] Morrison | Design and testing of a multi-sensor pedestrian location and navigation platform[END_REF], developed at the University of Calgary. The platform includes a Novatel receiver and supports up to ten 6DoF Analog Devices ADIS16375 IMUs. All data were synchronized with GPS time. During these field tests, two IMUs were used: one placed in the hand and one mounted on the user's foot, as shown in the right side of Figure 7. The foot mounted IMU served as a reference for the assessment of the step detection process. Five women and five men were equipped with this hardware setup and requested to walk twice along a curved route of about 300 m, for a total of 600 m. During the first run, the subjects were holding the inertial sensor in a swinging hand and for the second run, they were asked to change to the texting mode without stopping their walk. This corresponds to a natural change of the sensor carrying mode without stopping the walk. Because the pedestrian's walk is different for each person, the set of parameters K can be optimized individually through a calibration phase. To assess the benefit of calibrating the model, each subject was also requested to perform two straight paths in texting and swinging modes. This part of the test was used to find a calibrated set of K for each subject. For this purpose, -true‖ step lengths were evaluated by interpolating post-processed differential GPS positions over each identified step. Finally iterative least-squares (LS) process was applied to determine the calibrated set K n , with a convergence achieved at the n th iteration.

Experimental Results

In order to assess the performance of the proposed step length model, the user's motion mode and step events are first identified. The performances of the motion mode classifier, which are evaluated for the soccer field tests, are reported in Table 1. The rows of the table show the tested activities, namely swinging and texting modes, while the columns indicate the predicted motion modes. The confusion matrix reports along its principal matrix the percentages of correct detection for each state while the off diagonals report the percentages of misclassification. In Table 2 the probability of correct detection of the user's motion mode is also reported for each test subject along with the percentage of correct step detection. However an extensive validation of the classifier and of the step detection algorithms has been proposed in [START_REF] Susi | Motion mode recognition and step detection algorithms for mobile phone user[END_REF]. Then, the percentage of error over the travelled distance has been computed for each test subject with:

handheld GPS d d   (12)
where d handheld is the total traveled distance estimated using the step length model and d GPS is the reference distance evaluated using post-processed GPS carrier phase signals in a differential mode. The accuracy of the post-processed solution with a 1 km baseline in the open sky was better than 1 cm. A pedestrian dead reckoning approach has been used to compute the total traveled distance. It is based on the sum of each step displacement vector p t , estimated at the time t and starting from the initial position p 0 as:

 

where s t is the step length evaluated by applying the proposed model and θ t is the GPS-based walking direction over one step. This angle has been extracted from the GPS trajectory, post-processed in differential mode, for translating the estimated displacement information into the positioning domain.

As seen mathematically in Equations ( 12) and ( 13), using these post-processed headings does not affect the estimation of the distance travelled error. The error percentages are reported in Table 1 for all test subjects using both the universal and the calibrated set of parameters defined in Equation [START_REF] Kim | A step, stride and heading determination for the pedestrian navigation system[END_REF].

Men are indicated with -M‖ and women with -W‖. The number of iterations necessary to achieve the LS convergence in the calibration process is also reported. Results show that even for the universal model, most of the travelled distance errors are between 4 and 6%. Two larger percentages are obtained for the male subjects M4 and M5 with 8 and 9% respectively for the universal model. However, with the calibrated handheld based step length model, the results significantly improve and the highest percentage becomes 5%.

Overall, the performance of the handheld based step length model is comparable to those of models proposed for body fixed sensors [START_REF] Jahn | Comparison and Evaluation of Acceleration Based Step Length Estimators for Handheld Devices[END_REF][START_REF] Alvarez | Comparison of Step Length Estimators from Wearable Accelerometer Devices[END_REF]. The quality of the universal model is further confirmed by the low number of iterations required to reach convergence. In order to further assess the validity of the universal model, the absolute differences between fitted and universal parameters have been evaluated for each test subject. In Figure 9, minimum, maximum and mean values of the above differences show the variations between fitted and universal parameters. It is worth mentioning that the test subject M5 was significantly taller than most of the subjects that participated in the training phase. This explains the larger disagreement between the model and the truth. Increasing the number and the variety of training subjects would improve the performance of the universal model.

Conclusions

In this paper, a step length model for evaluating the distance travelled by a pedestrian holding an IMU in a hand has been proposed. Before computing the distance walked over a step, the carrying mode of the device, i.e., texting, phoning or arm swing, and the user motion are identified. By using adaptive algorithms, step events are detected from the signals recorded in the hand and for all arm/hand motions. The proposed step length model for smart phone users combines the user's step frequency, the user's height and a set of three variables. The algorithms developed for estimating the pedestrian's step frequency using non-body fixed sensors have been presented. The fact that the strongest frequency of the signal extracted from the handheld IMU can be coupled either with the step or the stride of the pedestrian is used to estimate the step cadence. Performance of the model was assessed in the position domain with a universal set of parameters and one fitted for each person. The experimental tests have demonstrated percentages of error over the travelled distance between 2.5% and 5%. The latter are similar to those achieved in the literature but with sensors rigidly attached to the body. These results are applicable to autonomous navigation and tracking of pedestrians using smart phones. In addition, the proposed algorithms are processing the IMU signals in windows between 1.5 to 2.5 seconds, which enables real time implementation. Further improvements of the step model could be obtained by increasing the typology of subjects involved in the model's training phase.

  pre-processed components of the accelerometer and gyroscope signals, respectively. This feature is used to distinguish the static from the dynamic state.

Figure 1 .

 1 Figure 1. Spectrogram of the accelerometer signal for the sensor mounted on the walking user's foot.

Figure 2 .

 2 Figure 2. Spectrogram of the accelerometer signal for the sensor placed in the walking user's swinging hand.

Figure 3 .

 3 Figure 3. Block diagram of the algorithm for the user's step length estimation using a handheld sensor.

Figure 4 .

 4 Figure 4. Normalized PSD of the accelerations sensed by the foot mounted sensor and the one in the swinging and texting hand. For both motion modes, the strongest frequency is coupled with the step event.

Figure 5 .

 5 Figure 5. Normalized PSD for the accelerometer sensed by the user's swinging hand.Here, the strongest frequency is coupled with stride events.

Figure 6 .

 6 Figure 6. Estimated, true step lengths and step frequencies computed with signals from a handheld IMU when the user is walking with his hand swinging.

Figure 7 .

 7 Figure 7. (Left side) Data collection set up for training the parameters of the step model. The subject walks at different speeds with one IMU in the hand and one on the foot. A second person paces the test subject by using a wheel speed sensor; (Right side) Data collection set up for assessing the step length model. The subject walks at his preferred speed with two (only one is used for data analysis) IMUs in the hands and one on the foot. The NavCube is carried at the waist.

Figure 8 .

 8 Figure 8. Linear fitting of the true step lengths (blue dots) with the user's heights (h) and the product of the strongest dominant frequencies with the user's heights (hf) at different walking speeds and hand's motions. The outcome is the universal set K in the step length model Equation (8).

Figure 9 .

 9 Figure 9. Minimum, mean and maximum absolute differences between -fitted‖ and -universal‖ parameters of the proposed step length model.

Figure 10 .

 10 Figure 10. Reference path in green and estimated trajectories: modeled step length with the universal parameters in red and with the calibrated parameters in blue for the test subject with the worst performance (M5).

Table 1 .

 1 Confusion Matrix for the motion mode classifier.

	Test Motions	Texting	Classified as: Swinging	Irregular
	Texting	100%	0%	0%
	Swinging	0%	98%	2%

Table 2 .

 2 Metrics for Evaluating the Handheld Based Step Length Model.

	Subject	%P det (motion)	%P det (steps)	%DistanceTravelled Universal Fitted Model Model	Convergence Iterations
	M1	100	99	5.8	5	4
	M2	100	100	4.8	4.3	3
	M3	99	100	5	4.5	3
	M4	100	99	8	4.2	6
	M5	100	100	9	3.8	7
	W1	98	97	5.2	4.3	4
	W2	100	100	3.2	2.5	3
	W3	98	99	4.5	4	3
	W4	100	98	5.6	5	3
	W5	100	100	5.8	5	4
	Mean	99.6	99.2	5.7	4.2	4

Acknowledgments

The financial support of Research In Motion (RIM), the Natural Science and Engineering Research Council of Canada, Alberta Advanced Education and Technology and Western Economic Diversification Canada is acknowledged. The authors would like also thank Aiden Morris for the design and the assembly of the NavCube and for his support during the test fields.


Sergio Lescano
email: sergio.lescano@femto-st.fr

Dimiter Zlatanov
email: zlatanov@dimec.unige.it

Micky Rakotondrabe

Nicolas Andreff
email: nicolas.andreff@femto-st.fr

Kinematic Analysis of a Meso-Scale Parallel Robot for Laser Phonomicrosurgery

Keywords: Kinematic, Parallel Kinematic Machine, Screw Theory, Piezoelectric, Phonosurgery

The paper presents the kinematic model of a new meso-scale (∼1 cm 3 ) parallel kinematic machine intended for laser phono-microsurgery of the vocal folds. The proposed mechanism uses the displacement generated by piezoelectric cantilever actuators manufactured via a Smart Composite Microstructure technique. The architecture, the geometry, and the position kinematics of the device, modeled as a spatial linkage, are discussed briefly. Then, the paper presents a velocity and singularity analysis and concludes that the new meso-scale parallel kinematic machine does not have singularities in the neighborhood of the reference configuration where it is required to operate.

Introduction

Laser phonosurgery is a surgical procedure used to enhance the voice using a laser beam onto the vocal fold. This paper describes a meso-scale parallel kinematic machine (PKM) conceived with the aim to improve the quality of laser phono-surgery. This mechanism is placed at the end of a laryngoscope and orients a mirror used to direct the laser beam onto the vocal fold. In order to account for all scale effects, to satisfy the constraints of the microworld, and to ensure the required accuracy, piezoelectric material [START_REF] Ballas | Piezoelectric Multilayer Beam Bending Actuators: Static and Dynamic Behavior and Aspects of Sensor Integration[END_REF][START_REF] Clévy | Signal Measurement and Estimation Techniques for Micro and Nanotechnology[END_REF] combined with a smart composite microstructure (SCM) fabrication process [START_REF] Wood | Microrobotics using composite materials: the micromechanical flying insect thorax[END_REF][START_REF] Wood | Microrobot design using fiber reinforced composites[END_REF] are proposed to build this meso scale PKM.

In Section 2, the kinematic requirements of the meso-structure are described, the mechanism's architecture and geometry are proposed, and the position kinematics is discussed in brief. Section 3 is devoted to velocity and singularity analysis using screw theory. The main conclusions are drawn and possibilities for future work are outlined in Section 4.

2 Kinematic Requirements and Proposed Structure

Kinematic Requirements

Adjus table Laryn gosco pe

The PKM is located at the head of an adjustable laryngoscope. The laryngoscope's diameter is 20 mm, and the space assigned for the robot is 10×10×10 mm 3 . The end of the laryngoscope is introduced into the cavity of the human larynx, at a distance of up to 20 mm from the vocal fold, as depicted in Fig. 1. To guide the direction of the laser beam onto the vocal fold, the PKM needs to have two rotational degrees of freedom (2 DOF). The standard average length of an adult vocal fold is between 17 and 20 mm. With a working distance of up to 20 mm, the PKM can scan the entire vocal fold if it can control the laser direction in a cone with a half-angle of 15 • . The SCM manufacturing technique mentioned above, allows us to build the meso-scale robots as was shown by [START_REF] Wood | Microrobotics using composite materials: the micromechanical flying insect thorax[END_REF][START_REF] Wood | Microrobot design using fiber reinforced composites[END_REF], first in a plane and then folding and popping out to form the 3D structure.

Minimal bandwidth of 200 Hz and a scanning resolution better than 100 µm are required to prevent overheating and damaging healthy tissue. Despite their small displacement, piezoelectric cantilever actuators have been chosen for their nanoresolution, work speed, and adaptability to the required scale [START_REF] Ballas | Piezoelectric Multilayer Beam Bending Actuators: Static and Dynamic Behavior and Aspects of Sensor Integration[END_REF][START_REF] Clévy | Signal Measurement and Estimation Techniques for Micro and Nanotechnology[END_REF]. Piezoelectric cantilever actuators with dimensions of 10 mm × 2 mm reach displacements of 100 µm. Piezoelectric cantilevers are arranged in parallel to not surpass the assigned dimensions.

Proposed Mechanism Architecture and Geometry

The proposed PKM structure is depicted in Fig 2 . The end-effector, or the moving platform, is used to orient a mirror, Fig 1.

The platform is connected to the base (attached to the laryngoscope) by a passive U-joint at P (leg 3), which constrains it to only tilt and pan. We define the base reference frame Pxyz with coordinate unit vectors i, j, k, and a rotating frame, fixed to the platform, with origin at P and coordinate unit vectors u, v, w, Fig. 2.

The 2-dof platform motion is actuated with two identical 6-dof RUS legs, legs 1 and 2. Piezoelectric actuation is used and modeled as the first R-joint of the each actuated leg, i.e., the piezoelectric cantilever assumed to be a rigid rod rotating about a fixed axis. The two active joints have a common axis, A 2 A 1 , with unit vector

k 1 = √ 2/2i - √ 2/2j. (Point A i is the projection of the U-joint center B i on this axis.)
In each active leg, the axis of joint 2, the first joint of the passive U-joint centered at B i , i = 1, 2, is parallel to A 1 A 2 , while the center, C i , of the spherical joint is along the common perpendicular of the U-joint axes.

At the reference configuration, when the angles of the U-joint at P are

θ 1 = θ 2 = 0, we have c 1 = --→ PC 1 parallel to u = i and c 2 = --→ PC 2 parallel to v = j.
Moreover, at the zero configuration, the distal links, B i C i , are parallel to the z axis, while the proximal links, A i B i , are perpendicular to it and at equal angles with the x and y axes.

Each joint of the PKM is a compliant hinge and is made of a polyimide film which allows great flexibility without snapping [START_REF] Lobontiu | Compliant Mechanisms: Design of Flexure Hinges[END_REF], Fig 3.

The dimensions of the PKM are given by |A i B i | = 20l (corresponding to the piezoelectric-cantilever actuator), |B i C i | = 2l, |PC i | = l, where l = 0.5 mm. The distance between Pz and the common actuated-joint axis,

A 1 A 2 , is (20 - √ 2
2 )l, while points A 1 and A 2 are at √ 2l from each other. Therefore,

-→ PA 1 = -(10 √ 2 -1)li -10 √ 2lj -2lk -→ PA 2 = -10 √ 2li -(10 √ 2 -1)lj -2lk (1)

Position kinematics

Given a feasible platform orientation, defined by the values of the passive-leg Ujoint angles, (θ 1 , θ 2 ), the mechanism configuration can be calculated by finding the intersection of a sphere centered at the point C i , i = 1, 2, with radius |B i C i | and the circle with center A i , with radius

|A i B i |, in a plane normal to line A 1 A 2 .
In general, there will be four solutions of the inverse kinematics. Points C i are obtained as:

--→ PC 1 = lu = c 2 li + s 1 s 2 lj -c 1 s 2 lk (2) --→ PC 2 = lv = c 1 lj + s 1 lk (3)
where c 1 = cos θ 1 , s 1 = sin θ 1 . Indeed, the platform orientation is given by

R(θ 1 , θ 2 ) = R x (θ 1 )R y (θ 2 ) =   c 2 0 s 2 s 1 s 2 c 1 -s 1 s 2 -c 1 s 2 s 1 c 1 c 2   (4)
We have

--→ B j C j = --→ A j C j - --→ A j B j
, and so

| --→ B j C j | 2 = 4l 2 = | --→ A j C j | 2 + 400l 2 -2 --→ A j C j • --→ A j B j (5)
where

--→ A j B j = 10l( √ 2c j 1 i + √ 2c j 1 j + 2s j 1 k)
, and c j i = cos q j i , s j i = sin q j i . Using ( 1)-( 3) and ( 5) it is now possible to obtain a linear equation in terms of c j 1 , s j 1

F j c j 1 + G j s j 1 = H j (6)
where F j , G j , and H j depend only on θ 1 , θ 2 , resulting in two solutions for q 1 1 . E.g. for θ 1 = 15 • , θ 2 = -15 • the solutions for (q 1 1 , q 2 1 ) are:

(0.720102952001288 • , 0.742258398498497 • ) (0.720102952001288 • , 12.160623836841685 • ) (12.16311819616120381 • , 0.742258398498497 • ) (12.16311819616120381 • , 12.160623836841685 • ) (7)

Velocity and Singularity Analysis

The velocity analysis is performed using screw theory, which has found a wide application, in particular in the the study of parallel robots. Instantaneous motion and forces acting on a rigid body can be represented by a twist and a wrench, respectively. With respect to a given reference frame, a twist (or a wrench) is given by a pair of (non-commensurate) physical vectors, (ω ω ω, v), the body's angular velocity and the linear velocity at the origin (or (f, m), the resultant force and moment with respect to the origin). A normalized twists or wrenches has the same 3D geometrical representation, called a screw, defined by an axis and a pitch. (For a translation or a couple, there is no axis, but only a direction, and the pitch is said to be infinite.)

The reciprocal product, ω ω ω • m + v • f, measures the power exerted by a wrench on a twist. Two screws, and any wrench-twist pair on them, are said to be reciprocal to each other if the product is zero.

Velocity kinematics

On a PKM the instantaneous twist, ξ ξ ξ , of the platform with respect to the ground can be expressed as a linear combination of the joint screws of leg j [START_REF] Mohamed | A Direct Determination of the Instantaneous Kinematics of Fully Parallel Robot Manipulators[END_REF]. Denoting by q j i and ρ ρ ρ j i the amplitude and the unit zero-pitch twist joint twist associated with the i th joint of the j th leg, we have:

ξ ξ ξ = q1 1 ρ ρ ρ 1 1 + q1 2 ρ ρ ρ 1 2 + q1 3 ρ ρ ρ 1 3 + q1 4 ρ ρ ρ 1 4 + q1 5 ρ ρ ρ 1 5 + q1 6 ρ ρ ρ 1 6 (8a) ξ ξ ξ = q2 1 ρ ρ ρ 2 1 + q2 2 ρ ρ ρ 2 2 + q2 3 ρ ρ ρ 2 3 + q2 4 ρ ρ ρ 2 4 + q2 5 ρ ρ ρ 2 5 + q2 6 ρ ρ ρ 2 6 (8b) ξ ξ ξ = θ1 ρ ρ ρ 3 1 + θ2 ρ ρ ρ 3 2 (8c)
Since legs 1 and 2 have 6 dof they impose no constraint on the platform. Therefore the possible twists, ξ ξ ξ , at every configuration, are given by (8c) as any linear combination of the passive leg's two revolute-joint unit twists, ρ ρ ρ 3 1 and ρ ρ ρ 3 2 . When such a well-defined continuously-changing basis exists for the twist system of the platform, the velocity kinematics can be expressed with as few equations as the parallel mechanism's dof [START_REF] Zoppi | On the velocity analysis of interconnected chains mechanisms[END_REF].

For this, one can first substitute (8c) into (8a) and (8b), and then eliminate the passive joint velocities using reciprocal screws obtaining a system of two equations relating the platform angular speeds θ1 and θ2 with the input velocities q1

1 and q2 1 . For each active leg, there is a unit pure force, ϕ ϕ ϕ j , j = 1, 2, reciprocal to all five passive-joint screws. Its axis passes through the centers of the leg's U-and S-joints.

(In a leg-singular configuration, where the center of the spherical joint is in the plane of the U-joint, there will be a whole pencil of concurrent planar forces with this property, and two, rather than one, reciprocal screws for that leg must be used to obtain the input-output velocity equations. We will assume that such leg postures do not occur during the operation of the mechanism.)

Taking the reciprocal product of ϕ ϕ ϕ 1 with (8a) and of ϕ ϕ ϕ 2 with (8b), we obtain:

ϕ ϕ ϕ 1 • ρ ρ ρ 3 1 ϕ ϕ ϕ 1 • ρ ρ ρ 3 2 ϕ ϕ ϕ 2 • ρ ρ ρ 3 1 0 θ1 θ2 = ϕ ϕ ϕ 1 • ρ ρ ρ 1 1 0 0 ϕ ϕ ϕ 2 • ρ ρ ρ 3 2 q1 1 q2 1 (9)
In ( 9), the term ϕ ϕ ϕ 2 • ρ ρ ρ 3 2 is zero, because the joint-2 axis of the passive leg always passes through point C 2 . To calculate the other entries, we can use:

ρ ρ ρ 3 1 = i 0 , ρ ρ ρ 3 2 = v 0 , ϕ ϕ ϕ j = b j |b j | c j × b j |b j | , ρ ρ ρ j 1 = k 1 -→ PA 1 × k 1 , j = 1, 2 (10)
I Since -→ PA j = c j -b j -a j we have that ϕ ϕ ϕ j • ρ ρ ρ j 1 = b j |b j | • (c j -b j -a j ) × k j 1 + k j 1 • (c j × b j |b j | ) = b j |b j | • [k j 1 × a j ] j = 1, 2 (11)
and so (9) can be written as:

  b 1 |b 1 | • k 1 × a 1 0 0 b 2 |b 2 | • k 1 × a 2     q1 1 q2 1   =   c 1 × b 1 |b 1 | • i c 1 × b 1 |b 1 | • v c 2 × b 2 |b 2 | • i 0     θ1 θ2   (12)
which can be abbreviated as B q = D θ θ θ , where B is a diagonal matrix.

Singularity Analysis

Gosselin and Angeles suggested a classification of singularities for parallel manipulators into three main groups associated with the singularity of square matrices such as D and B above [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF]. Later Zlatanov gave a more comprehensive study of the singularities of closed-loop mechanisms including configurations where the matrices D and B are not square [START_REF] Zlatanov | Analysis of the instantaneous kinematics and singular configurations of hybrid-chain manipulators[END_REF][START_REF] Zlatanov | Identification and classification of the singular configurations of mechanisms[END_REF]. This can occur in the case mentioned above of either leg 1 or leg 2 being in a singular posture with the center of the S-joint lying in the plane of the U-joint, in which case a redundant passive mobility singularity is present. However, if we assume that such leg postures are outside of the range of operation of the mechanism, the singularities of the mechanism are the configurations where either B or D is singular, leading to a degeneration of the inverse or the direct instantaneous kinematics, respectively. In ( 9) or ( 12), matrix B is singular when one or both elements of its diagonal becomes zero, signifying that the actuated joint axis intersects one of the reciprocal force axes. In this case, if D is nonsingular, we have a singularity of type (IO, RI), where the end-effector loses freedoms and the input velocities can be nonzero with a stationary platform.

This cannot happen in a neighborhood of the reference configuration where q j 1 = 0, matrix B is equal to diag(|a 1 1 |, |a 1 2 |) far from having a zero determinant. To find a singularity of the matrix D, the axis of either ϕ ϕ ϕ 1 or ϕ ϕ ϕ 2 must lie in the plane of the passive-leg U-joint. In this case, the platforms freedoms cannot be controlled. It will be able to move with actuators locked and the actuated velocities are not independent and cannot be chosen arbitrarily, a singularity of types RO and II. Here again, with a good choice of geometric parameters, the PKM is safe in a neighborhood of the position where q j 1 = 0 with

D = 0 -|c 1 | |c 2 | 0 (13)
Since the PKM workspace is limited to the angles required to scan entire vocal fold (±15 • ), both B and D are expected nonsingular, in this reduced workspace.

Conclusion

This paper presents a meso-scale PKM dedicated to the laser phonomicrosurgery of vocal folds. A kinematic analysis is performed. Reciprocal screws are used to find the two matrices that describe the mechanism's velocity behavior. All possible cases of singularities are analyzed. The limitation of the workspace and the proposed mechanism geometry avoid each type of singularity.

In the future, owing to the microfabrication limits and the robot structure compliance, a study of the directional variation of the screws will be taken into account using interval analysis [START_REF] Merlet | Interval analysis and reliability in robotics[END_REF]. Moreover, due to the low displacement generated by piezoelectric cantilever actuators, a different configuration in a position near-singular for the matrix D could be exploited [START_REF] Voglewede | Two Viewpoints on the Unconstrained Motion of Parallel Manipulators at or Near Singular Configurations[END_REF] to generate large platform rotations with small displacements of the actuators, being near the point of singularity without reaching it to preserve controllability.

Fig. 1

 1 Fig. 1 PKM, at the head of the laryngoscope, driving the laser beam direction onto the vocal fold.

Fig. 2

 2 Fig. 2 Platform operated in parallel by two active legs and constrained by an unactuated U-joint leg.

Fig. 3

 3 Fig. 3 Realization of the passive U-joint.

Acknowledgments

This work has been supported by the European Union Seventh Framework Programme [FP7- ICT-2011-7]. Project µRALP (www.microralp.eu).


Jean-Philippe Tarel

Rachid Belaroussi

ROBUST 2D LOCATION OF INTEREST POINTS BY ACCUMULATION

Keywords: Interest point, corner, feature point, detection, extraction, accumulation transform, repeatability, edges

Various interest point and corner definitions were proposed in the past with associated detection algorithms. We propose an intuitive and novel detection algorithm for finding the location of such features in an image. The detection is based on the edges in the original image. Interest points are detected as accumulation points where several edge tangent lines in a neighborhood are crossing. Edge connectivity is not used and thus detected interest points are robust to partial edges, outliers and edge extraction failures at junctions. One advantage of the approach is that detected interest points are not shifted in location when the original image is smoothed compared with other approaches. Experiments performed on Oxford and Cambridge reference databases allow us to show that the proposed detection algorithm performs better than 9 existing interest point detectors in terms of repeatability from multiple camera views.

INTRODUCTION

Interest points and interest point neighborhoods are highly used features to summarize an image into a very reduced set of information. For instance, interest points are of great importance for object tracking, object recognition, camera calibration and image indexing.

Given an image of a complex scene, it is hard for people to tell where the corners are. The selected positions differ from one person to another. This explains why several corner definitions, applying to images, were proposed in the past. Due to the difficulty to agree on reference corner positions, corner detection algorithms were difficult to rate [START_REF] Mokhtarian | Robust image corner detection through curvature scale space[END_REF]. Evaluation based on the repeatability of the detected points in different views [START_REF] Schmid | Evaluation of interest point detectors[END_REF][START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF] led to the introduction of interest points and brought reliable results.

Two classes of interest point detectors can be identified: pixel-and edge-based approaches. A recent example of pixelbased algorithm is [START_REF] Rosten | Faster and better: A machine learning approach to corner detection[END_REF], and an edge-based algorithm is de- scribed in [START_REF] Awrangjeb | A fast corner detector based on the chord-to-point distance accumulation technique[END_REF]. The here proposed algorithm is based on accumulation and thus belongs to the second class. However, contrary to most edge-based algorithms, it does not rely on edge connectivity, so it is more robust to missing edgels (edge elements), outliers and edge extraction failures at junctions. Another advantage lies in its ability to maintain accuracy location even when the original image is smoothed. The proposed algorithm can be used in many tasks, for instance for road sign detection [START_REF] Belaroussi | Angle vertex and bisector geometric model for triangular road sign detection[END_REF].

After the detailed description of the interest point detection algorithm based on accumulation, a comparative evaluation with state of the art interest point detectors is given, showing the benefit of the proposed algorithm.

DETECTION BY ACCUMULATION

Definition of interest points location

Our idea is to define interest points as points where several tangent lines to the edges are crossing. Our definition is based on the image edges and thus requires edge extraction as a first step. But contrary to many existing methods based on edges, such as [START_REF] Awrangjeb | A fast corner detector based on the chord-to-point distance accumulation technique[END_REF], we do not want to use the edge connectivity property. Indeed, edge connectivity is often compromised by various perturbations such as noise, occlusion, outliers and smoothing. The second step is an accumulation step over selected pairs of edgels. Third step allows finding the interest point locations and eventually sorting the obtained interest points in decreasing order of strength.

With the previous definition, an interest point is detected at its correct position, even when the corner is smoothed or truncated. This is illustrated at the top and bottom right corners of the images in Fig. 1. It is also important to notice that when intensity smoothing is applied on the image, a detected interest point is not affected by a shift in position like in many other approaches. These interesting properties of robustness and accuracy are a consequence of the accumulation of evidence on the crossing of edge tangent lines. In addition to the missing edgels robustness, the process is also robust to outliers as illustrated in the top left corner of the images in Fig. 1. Notice that points not usually considered as corners are also detected as interest points. For instance, in the bottom left interest point of the images in Fig. 1, two interest points are detected as the crossing of horizontal and vertical edges even though the edges do not actually intersect.

Edge extraction

Whatever the chosen edge extraction algorithm, it must provide a set of edgels. Here, an edgel is a 2D point position P i in pixels with the gradient vector G i perpendicular to the edge at P i . Usually, the edge extraction algorithm is parameterized by the kernel size s of the filter used to reduce the noise on image intensities. Edge detection is applied on the original image after normalization between [0, 255]. In practice, we use Deriche edge detector [START_REF] Deriche | Using Canny's criteria to derive a recursively implemented optimal edge detector[END_REF] which gives the horizontal and vertical components of the gradient at each pixel position in the original image I(P ). Then, a no-maximum suppression is applied to obtain edges. Edgels with a gradient norm lower than a threshold g m are discarded. This threshold allows removing edges in noisy areas and reducing the number of edgel pairs used in the next step. The output of the first step is thus a set of n edgels (P i , G i ).

Accumulation over edgel pairs

To find the points where several tangent lines are crossing, pairs of edgels with indexes i and j are selected. To reduce the set of possibilities and because interest points are defined locally, only pairs of edgels at a distance d ij lower than d m are selected. Moreover, to avoid detecting interest points along a low curvature edge, a selection is performed with respect to the angle between gradient vector G i and the opposite of gradient vector G j . We denote this angle by α ij , see Fig. 2, and it must be within [-π 2 , π 2 ] for selected edgel pairs. As a consequence, interest points with obtuse angle will not be detected. To take into account the estimation error on gradient angles, the interval

[-π 2 , π 2 ] becomes [-π 2 -α m , π 2 + α m ]
where α m is the assumed maximum error on gradient angles.

For a pair of selected edgels (P i , G i ) and (P j , G j ), we compute, when it exists, the crossing point C ij of the two tangent lines of the edges going through P i and P j respectively, see Fig. 2. From the equation of the tangent lines, we deduce two equations on the position C ij : G t i (C ij -P i ) = 0 and G t j (C ij -P j ) = 0, where t denotes the transpose operator.

C i,j

P j d i,j P i G i α i,j
G j Fig. 2. Accumulation at the crossing C i,j of the two tangent lines going through P i and P j with normals G i and G j respectively. The crossing point C ij is thus obtained as:

C ij = (G i G j ) -t G t i P i G t j P j (1)
The positions of C i,j over the selected pairs of edgels are accumulated in a 2D histogram h(P ) the size of which is the same as that of the original image I(P ). Each vote at position C i,j can be simply an increment by one. However, to favor edges with strong gradients over edges with low gradients, we use a vote with the value G i G j . This value is the logarithm average of the gradient norms G i and G j . In summary, the 2D histogram h(P ) can be formally described as:

h(P ) = n i=1 n j=1 G i G j 1I Pi-Pj <dm 1I |αi,j |< π 2 +αm 1I P =Cij (2)
where 1I is the indicator function (if x is true 1I x = 1 else 0), α i,j = G i , G jπ, and C ij is given by (1).

Local maxima

After the 2D histogram h(P ) is obtained, the positions of local maxima give the interest point locations. A local maximum in the 2D histogram h is a pixel position P where h(P ) is higher than or equal to the values at its eight neighbors. The detected interest points can be sorted using the values of h at the maximum. The N first sorted points can be selected as the detected interest points. An alternative to sorting is simply to remove the locations where the value of h is smaller than a threshold.

EXPERIMENTS

Parameters selection

The gradient norm threshold g m used during edge extraction is usually between 10 and 40 gray levels, i.e between 4% and 16% of the maximum intensity of the input image. The higher the value of g m , the faster the interest point detection algorithm. The scale parameter s used to filter the image during edge extraction is the main parameter of the detection algorithm by accumulation since it specifies the scale at which interest points are detected. In noisy or randomly textured images, a large value of s should be preferred. On the contrary, in clean detailed images, a small value of s should suffice. Usually s is between 0.3 and 3 pixels. The maximum distance d m used for edgel pairs selection can be used to change the ordering of the detected interest points. But the higher the value of d m , the higher the number of edgel pairs and thus the slower the interest point detection algorithm. Moreover, d m must not be too small to allow stable interest points detection first. In practice, d m is set between 2% and 8% of the input image size. The last parameter α m does not play an important role and is usually set to 0.2 rad. It is also introduced in the algorithm to allow, when it is necessary, extracting interest points with obtuse angles as well as interest points with acute angles. In the left and right images of Fig. 3, the extracted interest points are shown in red, left with α m = 0.2 rad and right with α m = 1.3 rad. With α m = 0.2, only acute interest points are detected; when with α m = 1.3 both acute and obtuse interest points are detected. The price to pay for the detection of obtuse interest points is an increasing number of low stability interest points along edges with small curvatures. The original image comes from1 , see [START_REF] Mokhtarian | Robust image corner detection through curvature scale space[END_REF].

Comparison with other interest point detectors

For the evaluation of the interest point detection algorithm by accumulation in comparison with state of the art interest point detectors, we follow exactly the evaluation scheme described in details in [START_REF] Rosten | Faster and better: A machine learning approach to corner detection[END_REF], similar in its principle to the one described in [START_REF] Schmid | Evaluation of interest point detectors[END_REF]. The evaluation is based on the ǫ-repeatability which is defined as the ratio of the number of interest points repeatedly detected in a pair of images (of the same scene) over the number of interest points appearing in both images. ǫ is the maximum error accepted on the interest point location (typically ǫ = 5 pixels). As in [START_REF] Rosten | Faster and better: A machine learning approach to corner detection[END_REF], two image databases with ground-truth mapping between views are used: the Cambridge database2 and the Oxford database 3 . The Oxford database consists in 8 sequences of 5 pairs with a homography mapping between two views, whereas the Cambridge database consists in 3 sequences of 210, 182 and 56 pairs with a 3D mapping between two views. For algorithm comparison, the repeatability is plotted versus the number of detected interest points for different algorithms. The algorithms parameters are selected by optimizing the results on the "bark" dataset, then the same values are applied to other datasets. Only the scale s of the smoothing is modified for the Cambridge database due to better registration between views compared with the Oxford database. The repeatability curves obtained in [START_REF] Rosten | Faster and better: A machine learning approach to corner detection[END_REF] are shown for Harris, Harris-Laplace, Shi-Tomasi, SUSAN, DoG (used in Sift), FAST-9 and FAST-ER algorithms, along with the curve obtained when interest points are randomly sampled in the image. We added the two blue curves which are for the Radial Symmetry Transform [START_REF] Loy | Fast radial symmetry for detecting points of interest[END_REF] (RST) and the proposed accumulation based interest point detection algorithm (Accum). The results obtained in Fig. 4 on the Cambridge database equal those of the best FAST-9 and FAST-ER algorithms which are based on learning on this database, with the advantage that the algorithm based on accumulation is particularly simple to implement and does not need a learning step.

Fig. 1 .

 1 Fig. 1. Extracted edges, accumulation of tangent line crossings, and detected interest points in red superimposed to the original image. Used parameters are s = 1, g m = 32, d m = 16 and α m = 0.2.

Fig. 3 .

 3 Fig. 3. Detected interest points in red superimposed to the original image when α m = 0.2, left and α m = 1.3, right. Interest points with obtuse angles are not detected in the left image and are detected in the right image. The other parameters are s = 1, g m = 32, d m = 16.

Fig. 4 .

 4 Fig. 4. Comparison of the repeatability curves of different interest point detectors on the 3 sequences of the Cambridge database 2 . Used parameters are s = 1 3 , g m = 32, d m = 16 and α m = 0.2.

http://www.ee.surrey.ac.uk/CVSSP/demos/corners/originals.html

http://svr-www.eng.cam.ac.uk/˜er258/work/datasets.html

http://www.robots.ox.ac.uk/˜vgg/data/data-aff.html

Acknowledgments

Deep thanks to Edward Rosten for providing us, from [4], the curves for 8 of the algorithms which appear in Fig. 4 and5.

* Thanks to the ANR (French National Research Agency) for funding, within the iTowns-MDCO project.

Compared with the other algorithms, ours scores are higher in 6 sequences over the 8 of the Oxford database, see Fig. 5. In the Leuven sequence, where gamma transformations are applied on the image intensities, our results are within the average range of the other results. This indicates a lack of robustness to gamma transformations due to the g m gradient threshold. This can be fixed by adding a preprocessing step enforcing normalization of the average of the logarithm of the image intensity. In the Wall sequence, where images are taken with various viewing angles, our results are similar to the best results despite a large zoom with a constant smoothing scale s.

Average processing time for interest point detection with sorting is 624 ms on the Oxford images (max size 1000×700) and 417 ms on the Cambridge images (768 × 576) using a 1.1GHz Intel Core 2 Duo PC. Thus, it is similar in speed to the Harris detector.

CONCLUSION

We described an interest point detection algorithm based on the accumulation of evidence on the position of the crossing of pairs of edge tangent lines. This algorithm uses edges but is robust to missing edgels and outliers. In the comparative experiments, our interest point detection algorithm shows enhanced results over 9 state of the art algorithms. We point out that the scale of the filter during the edge extraction is the important parameter related to the scale at which interest points are extracted. We thus plan to extend our algorithm by introducing scale self estimation.