Semi-automatic annotation of the UCU accents speech corpus

Annotation and labeling of speech tasks in large multitask speech corpora is a necessary part of preparing a corpus for distribution. We address three approaches to annotation and labeling: manual, semi automatic and automatic procedures for labeling the UCU Accent Project speech data, a multilingual multitask longitudinal speech corpus. Accuracy and minimal time investment are the priorities in assessing the efficacy of each procedure. While manual labeling based on aural and visual input should produce the most accurate results, this approach is error-prone because of its repetitive nature. A semi automatic event detection system requiring manual rejection of false alarms and location and labeling of misses provided the best results. A fully automatic system could not be applied to entire speech recordings because of the variety of tasks and genres. However, it could be used to annotate separate sentences within a specific task. Acoustic confidence measures can correctly detect sentences that do not match the text with an EER of 3.3%

Introduction
The UCU Accent Project is a medium sized longitudinal multilingual multitask speech corpus, currently being collected at University College Utrecht (UCU), in the Netherlands. The student community is international and English is the lingua franca, providing a rich source of multilingual and multiaccented native and non-native speech data.
Three student cohorts are being recorded over a six year period and the final corpus size will be upwards of 1000 recordings, producing over 600 hours of speech. Each recording contains up to 12 speech tasks, including read texts of varying difficulty, monologues and dialogue. Speakers speak in both English and their L1.
For this paper, we have taken three approaches to segmenting and labelling the speech tasks -automatic, semi-automatic and manual -with a view to finding optimal efficacy and accuracy.
Manual labeling of speech data is still the most common form of labeling despite a large body of research on automatic methods. As observed by the many authors writing about manual and automatic labeling, for example, (Gut and Bayerl, 2004), (Campbell and Sagisaka, 1992), (Cole et al., 1994), (Truong and Trouvain, 2012), (François et al., 2012), there is much variability in accuracy and inter and intra annotator agreement, affecting the reliability and validity of the labels. Annotators can be affected by different L1 backgrounds, different dialects or differences in perceptual acuity. Even where there is agreement on the item to be labeled, there is the problem of interpretation of the labeling strategy and of the application of coding (Cosi et al., 1991). Manual labeling of large or longitudinal corpora also lends itself to a changing team of annota-tors. Individual familiarity with the project, motivation and interest in the task, and the extent and type of previous training will also influence annotator judgement. While manual segmentation and labeling might be expected to produce accurate results in single instances, the tedious and repetitive nature of the task leaves it vulnerable to fatigue and simple human error.
With such a large number of speech recordings to be labeled for specific speech tasks, we look at how much of this work can be automated, and whether automatic procedures can provide sufficient accuracy while reducing time investment and vulnerability to error. Among the automatic systems we looked at are acoustic event detection, and task and utterance segmentation based on automatic speech recognition (ASR).

The UCU Accents Speech Corpus
The UCU Accent Project is a longitudinal corpus. Four consecutive cohorts of students are being recorded at five time points over the three years of their undergraduate study. The majority of the speakers have Dutch as their L1. Apart from Dutch speakers and native speakers of English, over 30 other native languages are represented.
Facilitators include faculty members from UCU and from department of Linguistics at the Utrecht University, masters students and undergraduate students. All facilitators use the speech data from the recordings for their own research and as such are invested in maintaining the quality of the recordings.
The speaker is recorded on eight different channels, via different microphones, while seated at a table in one of the faculty offices in the college. A more detailed description can be found in (Orr et al., 2011).

Speaker Tasks
The tasks are varied in difficulty and type, and include both native and non-native speech. The speaker is asked to perform between 9 and 12 speech tasks, including texts that are read aloud, monologues in both the native language and English (if English is not the native language), and a dialogue with the person facilitating the recordings. The texts are intended for a variety of research purposes, and include the production of sentences for use in intelligibility testing, sentences for use in prosody analysis, and texts containing shibboleths. The speaker tasks are listed in Table 1. Extract from the Rainbow Passage (Fairbanks, 1960) 3 Please Call Stella * 4 The Boy who Cried 'Wolf' (Deterding, 2006) 5 Balanced sentence sets: intelligibility testing † 6 5 sentences for investigating rhythm (White and Mattys, 2007) 7 Extract: Declaration of Human Rights (L1) † † 8 Extract: Declaration of Human Rights (English) † † 9 2 minute monologue L1 10 2 minute monologue English (formal topic) 11 2 minute monologue English (informal topic) 12 3 minute dialogue with the facilitator * Speech Accent Archive, George Mason Univeristy http://accent.gmu.edu † Task 5 refers to sentences from van Wijngaarden (Van Wijngaarden et al., 2002) in quantifying intelligibility of speech in noise for non-native listeners. † † http://www.un.org/en/documents/udhr/

Variability of content in the speech files
The speech recordings contain quite some variability. Speech tasks were added as the corpus developed in order to make it comparable to other similar corpora. The order in which tasks were produced could also vary. Furthermore, the recordings were originally made with no clearly audible separation between tasks. After the first round of recordings, with concern for labeling the data, an audible separator was introduced between the tasks. The choice for an audible separator stems from the nature of the recording setup. It is not only a signal for separating tasks, but it is also a prompt for the participant, making clear when they should speak.

Semi-automatic acoustic event detection
The presence of an audible separator between tasks suggested the use of acoustic event detection to locate the beginning of each task and label it, producing a table of timestamps and corresponding events, i.e. speech tasks.
For this, we used a lightweight fully-connected twostate hidden Markov model with Viterbi decoding, with states for 'speech' and 'event,' using the implementation developed in (van Leeuwen, 2005). The output probabilities are modelled using 16-component Gaussian Mixture Models. These were trained on five manually segmented files from the recorded speech data, so that in total, approximately 50 tokens for 'event' were used for training, and approximately 40 minutes of speech. The transition probabilities were manually defined, in order to have maximum control over the false alarm/miss trade-off and class sequence smoothing. The transition probability matrix was parametrised by a parameter p trans and a prior over the states p bell which represents the prior that in Automatic Speech Recognition is modeled by the language model. Effectively, in decoding every frame the transition matrix is used to update the current maximum log likelihood path. We used a fixed parameter log p trans = −10 for smoothing detection transitions, and controlled the false alarm/miss trade-off via p bell . In manual checking of the detected events, it is much easier to dismiss a false alarm than to find a missing event. Hence we operated using a relatively high prior for the event.
As features, we used 13 MFCC 1 parameters plus deltas computed over 9 consecutive 25 ms frames taken at 10 ms intervals. We used the 'rastamat' implementation made available by Dan Ellis (Ellis, 2005) in GNU Octave (Eaton, 2002), using default HTK 2 compatible parameters. After feature extraction, an energy-based sound activity detection algorithm was employed, integrating energy in a 300-8000 Hz range. The criterion for frame selection was that the energy in this bandwidth was larger than 30 dB below the maximum energy in the recording. Non-sound frames were in effect removed from the modeling and the decoding, and later inserted in the time mark-up. This had the effect that sometimes the duration for the detected event appeared much longer than the actual event, because the Viterbi decoding would have spanned just over a larger chunk of silence. Since the detected events were only there to mark task boundaries, this was not considered a problem.

Evaluation
All detected events were manually checked. Events were played, and either accepted as instances of the 1 Mel Frequency Cepstral Coefficients 2 Hidden Markov Toolkit audible separator, or rejected as false alarm with a single key-stroke. If the number of accepted events was too low, a manual audio and visual scan through the sound file was made to find missed events.
For 292 recordings, the number of events was 3530 (a mean of 12.1 per recording), where 381 false alarms were generated and 13 events were missed. The total recording time excluding events for these recordings was 97.0 hours, with a total duration of events of 3.24 hours. Note that these times include silence that may arbitrarily be attributed to 'event' or 'speech'. Hence we were operating at 3.9 FA/hour and 0.37 % miss rate. Listening to and judging the events took 67 seconds per recording, on average, whereby the total "wasted time" on the false alarms was about 1100 seconds. This should be compared to the roughly 800 seconds necessary to manually locate the missed events in the audio recordings.

ASR-based task classification
For the first 115 recordings, where no audible task separator was used, we developed an ASR-based task classification system, inspired by the work in (Moreno et al., 1998). Without transcriptions for the spontaneous speech parts of the recordings, we cannot align a recognition of the entire recording on the transcription. In particular, because we use a language model trained solely on the transcriptions, fragments of the spontaneous speech task might be incorrectly aligned to the tasks in the transcription.
Therefore, we aligned part of the recognition to the transcriptions of each separate task. For each task with N words in the transcription we took the first 3 2 N words of the recognition and aligned the two texts. Then, we shifted the recognition with 1 4 N words and again aligned 3 2 N words with the transcription. We did this for each window of 3 2 N words until we reached the end of the recognition. In the end, the window with the highest number of aligned words was considered to be correct. We determined the exact start and end time of the task by selecting the first and last word that was aligned correctly. By using only marginally larger windows than the transcription length, we minimised the probability that parts the recognition of other tasks are aligned to the transcription.
We used our in-house developed Automatic Speech Recognition (ASR) system "SHoUT" (Huijbregts et al., 2009) with models for English and Dutch, without any further adaptation to the acoustic models. The entire recording was first pre-segmented in utterances based on energy and a minimum pause criterion, which worked well for the clean recordings in this corpus. Recognition results of all utterances in both languages were used as hypothesis transcription in the alignment procedure described above.

Evaluation
We have evaluated this automatic classification approach on a test set of 48 recordings. Table 2 contains for each task the percentage of time that the recording was falsely classified as the particular task (false alarms) and the percentage of time that the audio was not classified as the particular task when it should have been (missed speech).

Manual segmentation
The files without the audible separator were also segmented manually. The files were opened in a waveform editor 3 to allow playback for the listener as well as visual inspection. When a new task was located, a label was placed at a time point just before that task. The label file information was stored in the same format as the files generated by the event detection system.
Files were segmented in batches of 20 to avoid fatigue. In practise, however, because of the difference in the length of the read texts, and because of the clear differences in energy patterns between read texts and spontaneous speech, locating the speech task was very easy for all tasks with the exception of the switch from monologue to dialogue in some cases.
The manual task segmentation results can be used as a reference for further work with the automatic ASRbased approach.

Automatic utterance identification
The goal of the automatic utterance segmentation was to find the identity of all the intelligibility sentences of Task 5 automatically. Variability in the content of the speech files makes a strict force-alignment over the entire task impossible. Additionally, each sentence was scored on pronunciation accuracy to allow automatic or semi-automatic selection of the sentences that were pronounced correctly. We used forced recognition of entire sentences and a separate phone loop recognition to compare the acoustic likelihood of each recognised sentence with the likelihood in the free phone loop. We created a pronunciation dictionary in which the words in each sentence were concatenated to one single 'word' in the dictionary. Utterances were segmented using a simple silence detector, which works because the intelligibility sentences are short (8-9 syllables) and spoken in a single effort, with pauses between consecutive sentences. The language model in this task was a simple uniform unigram model. In this manner, each complete sentence can be recognised irrespective of order they are pronounced.

Evaluation
The result of the forced sentence recognition was a list of start/end times of the sentences with acoustic confidence and the identity of the sentence. In order to evaluate the accuracy of the sentence recognition, we manually checked the recognition of 2010 sentences in 48 recordings.
A detection-error trade-off plot (Martin et al., 1997) is shown in Figure 1. This shows the trade-off between false positives and false negatives in error detection. The equal error rate is 3.3 %, which is a single-valued summary of the detection performance. For the sample of 2010 sentences, 38 were actually pronounced wrongly. At 10 % false alarms, there are no more mispronounced sentences missed.
From the 2010 sentences in the test set, 98.96% was correctly classified as correct. 84.21% was incorrectly classified as correct.

Discussion and Conclusion
Of the different methods that we looked at for the main task segmentation task, a combination of event detection and manual checking was the most efficient in terms of accuracy and time investment. The average 1.1 minutes per speech file to manually check the results of the acoustic event detection compares favourably to the ca. 4 minutes that it took to label a file manually. The low-miss-rate operating point we have chosen led to comparable costs of misses and false alarms, being 800 and 1100 seconds, respectively. This usually is an indication of a good choice of operating point. Furthermore, the manual labeling could take quite considerably longer where the listener had to search for the change from monologue to dialogue, or to search in the balanced sentence sets.
The ASR-based task segmentation was not sufficiently accurate for our recordings. We have a clear preference now for working with the event detection system, alongside manual checking. A semi-automatic system seems to be the most efficient way of segmenting the tasks in this corpus.
The ASR-based utterance segmentation works well for the intelligibility sentence sets in Task 5. The manual checking script used in evaluation (cf. Section 3.4.) will be integrated into our labeling work flow, similarly to the event detection system.
The next steps include applying speaker diarization to the dialogue sections of the recordings, as well as further optimising the ASR-based task segmentation. Whether the ASR-based approach or the eventdetection approach proves more useful in the end, we will always have to run a manual check on the results before delivering the corpus for distribution.


A unified SLAM solution using partial 3D structure

1 Abstract —Good quality of environment mapping demands modelling the associated environment nearly to its 3D originality. This paper presents a unified Simultaneous Localisation And Mapping (SLAM) solution based on partial 3D structure. As compared to existing representations such as grid based mapping, the novelty of the proposed unified approach lies in estimation, representation and handling of compact partial 3D features-based map model for a team of robots that are working in an unknown environment with unknown poses. The approach replies on a camera to perceive the environment and a 2D laser sensor to generate a SLAM solution with partial 3D features based representation. Extended Kalman Filter (EKF) estimates the robot pose based on its motion model and map of the explored environment. The solution has been tested in an indoor environment on two identical custom-developed robots. Experimental results have demonstrated efficacy of the approach. The presented solution can be easily applied on a distributed/centralized robotic system with ease of data handling and reduced computational cost.

1 Abstract-Good quality of environment mapping demands modelling the associated environment nearly to its 3D originality. This paper presents a unified Simultaneous Localisation And Mapping (SLAM) solution based on partial 3D structure. As compared to existing representations such as grid based mapping, the novelty of the proposed unified approach lies in estimation, representation and handling of compact partial 3D features-based map model for a team of robots that are working in an unknown environment with unknown poses. The approach replies on a camera to perceive the environment and a 2D laser sensor to generate a SLAM solution with partial 3D features based representation. Extended Kalman Filter (EKF) estimates the robot pose based on its motion model and map of the explored environment. The solution has been tested in an indoor environment on two identical custom-developed robots. Experimental results have demonstrated efficacy of the approach. The presented solution can be easily applied on a distributed/centralized robotic system with ease of data handling and reduced computational cost.

I. INTRODUCTION
Simultaneous Localization and Mapping (SLAM) is considered as the fundamental problem in mobile robotics [1]. Initially both the map and the mobile robot position are not known and the robot moves through the environment, perceives it using on-board sensors and finally generates a map model. Due to sensor noises and process noises associated with motion of the robot, various probabilistic estimation techniques are used to explicitly model different sources of noises and their impact on measurements for efficient robotic map building.
The applications of SLAM extend from land environments to aerial and underwater scenarios [2]. A number of latest findings have been reported in the present decade. SLAM solutions have been successfully integrated with path planning and navigation issues to enhance autonomy of the moving robot. SLAM problem, at present, is considered as a solved problem with different estimation Manuscript received March 23, 2014; accepted July 18, 2014. This research was funded by grant from University of Genova, Italy.
approaches [1]. Many researchers have proposed distinct solutions starting from late 80's when stochastic techniques have started to replace deterministic based approaches in solving localization problem because of associated uncertainties in moving robots [3]. Due to variation in stochastic techniques, a number of SLAM versions have been found in the literature. Researchers have used Extended Kalman Filter (EKF), Information Filter (IF), Sparse Extended Information Filter (SEIF) and Rao-Blackwellised Particle Filter (RBPF) to propose SLAM solution [4]. Various combinations of integrated on-board sensors (e.g. 2D laser scanner, ultrasonic sensor, camera etc) have been used for testing such solutions [4]. In [5], EKF based SLAM solution has been presented using features for map modeling while SEIF based solution has been proposed in [6]. RBPF based approaches rely on map features modeling [7] or grids based map modeling [8].
In the present research, feature based EKF SLAM algorithm is proposed to solve localization and mapping issues. A 2D laser sensor and a camera are used to perceive the environment nearer to its 3D originality.
The organization of the paper is as follows: Section II discusses the related work. Section III highlights the proposed feature based representation scheme. Section IV presents unified EKF SLAM framework and its results have discussed in Section V. Finally Section VI comments on conclusion and potential applications of this research.

II. RELATED WORK
Scientific community reports several contributions to present 3D solutions using camera [9] since it is considered as cheapest and useful sensor to detect 3D features. Camera provides rich set of information about the environment which can be used for robot localization and mapping. Numerous techniques employed to use single camera or a group of cameras to efficiently solve SLAM problem are reported in [10]. Most of these approaches have used point features (using patches) extracted from the image. This approach has been considered useful; however it cannot efficiently produce a 3D model of the explored environment. Another major challenge offered by this approach is in the features initialization [9], [11]. Camera provides bearing A Unified SLAM Solution Using Partial 3D Structure S. Riaz un Nabi Jafri 1,2 , J. Iqbal 3 , H. Khan 1 , R. Chellali 1 1 only information of the feature and it is almost impossible to extract the depth information of the features using only one image (single view) [12]. An inverse depth concept has been used to initialize features with depth uncertainty reaching to infinity [12], [13]. It has been successfully demonstrated that this scheme works well not only to localize the camera but also to map very distant features. In an attempt to model lines using end points, use of line segments has been proposed in [14], [15]. Some drawbacks associated with this solution include delayed initialization and uncertain end points [9].
To solve SLAM problem, the present work exploits simplicity and computational efficiency of MonoSLAM i.e. single camera based SLAM solutions [16]. The distinguishing features of the proposed approach include building compact map model and sharing/ processing of information among team of robots.

III. PROPOSED APPROACH
In indoor environments, geometric shapes are common such as tables, walls, chairs and cupboards etc., which need to be modeled. In this work we are detecting these shapes from integrated sensors and extracting lines/corners features from them to model the environment.  horizontal Field of View (FoV) than the camera that is why we are observing some objects with horizontal features whose vertical counterparts are not detected. Though the camera has smaller horizontal view but its wide vertical sensing range permits it to observe window limits not detected by the laser. Figure 1(b) illustrates the desired representation of the environment with the robot (represented by triangle) and extracted features from both sensors. Laser sensor can display extracted feature immediately but it is almost impossible to display camera extracted feature exactly on the actual location because of unavailability of depth information. These limitations offer several challenges to use both sensors and demand a solution to initialize camera based features immediately and to converge them after getting sufficient parallax [12]. Here we are proposing a unified approach for presenting the features in a common way and to update them in EKF SLAM.

A. Camera Features
After image capturing from a calibrated camera, inverse depth parameterization of vertical lines has been incorporated in EKF framework. Consider a robot moving in an indoor environment from point A to point B with initial pose (xA ,yA , θA) to final pose (xB , yB , θB) as shown in Fig. 2. During its motion, it observes a vertical line feature (as shown by solid cyan line), which is projecting on X-Y plane at point V.
The bearing of the first vertical line feature is denoted as in the robotic frame of reference as shown by arrow from point A to point V in Fig. 2. The camera frame of reference and robotic frame of reference (shown by XR-YR axis with dashed arrows) are considered as same in this work. A 2D unit vector describing the direction of the line feature from the camera can be represented as (1) where is the u-axis (horizontal camera axis) central coordinate in the image, is the focal length of the camera, is the physical spacing of pixel on the optical sensor of the camera, is the rotation matrix from the camera coordinate frame to the world coordinate frame and is the detected u-axis coordinate of the vertical line feature. The bearing can be represented using directional vector as The world frame is shown by XYZ arrows. The bearing of the feature in world frame ( ) can be computed using The inverse depth representation for the projected feature point can be given by where ( , ) is the location of the robot in world frame, α is the bearing of the line in world frame and is the inverse of the range (depth) information of the line feature ( ) and is initialized with a guess 1/ .
Projected coordinates of V ( , ) can be written as During its motion from A to B, the robot is detecting again vertical line with different local bearing in its view as shown by arrow from B to V in Fig. 2. Point coordinates of V take form of (7) All the robot poses at A and B are known (with nearest global values) from the motion model. (6) and (7) show that only inverse depth ( ) is unknown and we need to find its value using significant changes of robotic pose as well as feature bearing angles. The proposed approach is based on taking the advantage of the robot successive movement and iteratively predicting a nearer estimate for the point V coordinates. For future matching in upcoming images, an image patch covering all vertical portion of the image with a small width around the line is stored to associate line in future observations. When the robot moves to B and takes an image then the line appears in the image on different horizontal pixel value. To associate this line with the older line, we need to predict a guess for the older line which shows that the older line will probably appear in this horizontal region of the new image shown in pixels. This guess is generated using motion model of the robot. Currently it is assumed that we are getting the pose prediction and = , which gives the location of the robot at point B in world frame. The older feature location can be predicted using (8) .
Rewriting (8) in inverse depth form as Equation (9) is used to find the horizontal pixel value of the older line for the image feature. Measurement model given by (10) shows the representation for finding the value of predicted older line position ( ) To associate new observed lines with the old ones, nearest neighbor method is used. For each new line, it is evaluated that new observation and old prediction satisfy (11) , p z u u u th   (11) where is the new observation for current feature and thu is the threshold (in pixels). After getting matched pair, for more robust association, new image patch belonging to a new current feature is compared with previously stored patch through cross correlation. In case at any step, if no line is matching with new line and new line is far from older lines (greater than predefined threshold) then new line is considered as fresh candidate and is initialized immediately.

B. Laser Features
This work uses laser based horizontal geometrical (line and corner) features and parameterize them as inverse depth features. Consider the same scenario of the robot movement as shown in Fig. 3. At point A, the robot observes a corner C (as pointed by brown arrow) with range (distance from robot) as and bearing as . The laser sensor frame of reference is fixed and is same as the robotic frame of reference. Bearing of the corner feature can be given by (12) .
The corners can be initialized and represented by   , , , , where symbolizes laser corner feature, ( , ) is the location of robot at point A, and are bearing (in world frame) and range of the corner from the robot respectively. Corner point coordinates C ( , ) in terms of (13) can be written as During motion from A to B, the robot is detecting again corner with different range and bearing (say = 5 ELEKTRONIKA IR ELEKTROTECHNIKA, ISSN 1392-1215, VOL. 20, NO. 9, 2014 ) as represented by arrow from B to C in Fig. 3. We need to predict old corner feature and to associate this with new observed corner feature. The prediction ℎ of old corner feature is given by where dhlc and αhlc are the predicted range and bearing for the old corner from point B respectively, (xB, yB, θB) is the predicted new robot pose at point B and (xC, yC) is the old corner coordinates as determined by (14). In this work, nearest neighbor method is used which is detailed in [17].
Each new corner measurement is tested with (old) prediction against threshold as given by . lc lc lc h z th   (16) If the re-observed corner is found within a predefined threshold then it firmly associates with the previously seen corner. If no match is found then the corner is treated as new feature and is initialized as given by (13).
At the starting location (point A), the robot is observing horizontal lines (in brown) shown in Fig. 4. One such line is pointed by brown arrow from point A in the figure. The detected line is making a global bearing of with range of . Feature initialization is given by where represents the laser line feature, ( , ) is the location of the robot and and are bearing and range (in world frame) of the observed line respectively. The end points global coordinates for the lines are saved to associate them in future. As the robot moves, the line feature will re observe with different range and bearing (say = ) as indicated by arrow from point B in Fig. 4. The prediction ℎ of old line feature is given by where dhll and αhll are the predicted range and bearing for the old line from point B respectively, (xB, yB, θB) is the predicted new robot pose at point B, and are the old line global range and bearing respectively. Each new line feature measurement is tested with old line feature prediction for validating threshold given by IV. UNIFIED EKF SLAM Based on fusion of features extracted from both kinds of sensors in the same algorithm, a unified EKF SLAM solution has been proposed.
The motion prediction only affects the robot state during motion without changing the features state. F and G are represent the Jacobians of motion function , , , w.r.t. system state and w.r.t. noise , , = , , ∆ respectively, Q is represents the covariance matrix for the noise , , as stated in [12].

C. Measurement Update
After prediction, in the measurement update step of the algorithm, Kalman gain K is determined by [3] where is Jacobian of the measurement prediction function w.r.t. the state vector. Three different functions for predicting different sensors as given by (10), (15) In (24), is the observation noise matrix, which deals with uncertainty of the observations [14]. After determining Kalman gain , state vector and covariance matrix are updated: where h and z represent the predictions associated with the features and current measurements respectively [12]. This procedure involves final state estimations and its covariance in the current iteration. New features are added after these steps. Though all kinds of features are distinguished w.r.t. different measurements and associated models but a balanced state and covariance updates result in a peaceful and unified solution.

V. EXPERIMENTS AND RESULTS
An indoor office like environment is selected for experiments having objects with orthogonal and nonorthogonal orientations. The robots are moved manually from different start-end positions to get different trajectories and corresponding laser range data and images have been recorded. The frequency of both sensor inputs is fixed to 10 Hz. Maximum linear speed of each robot is around 0.1 m/sec and maximum angular speed is around 0.03 rad/sec. The different data sets obtained are then used independently in EKF algorithm. The algorithm is developed by combining standard MATLAB functions with various open source functions such as openSLAM developed by research community. Hough transformation function is used to extract image lines while Split and Merge algorithm is used to extract laser lines.
In first experiment, the robot is initialized with pose (0, 0, 0). Figure 5(a) illustrates the first image from the camera recorded at the starting location with the vertical lines (shown in green) observed from the image. Figure 5(b) indicates the corresponding results of laser lines (in blue) and corner features (in cyan) extracted from the first scan. In this figure, vertical image lines representing inverse-depth projected points (in green with numbering) and current robot location (in red) are also shown. All the vertical lines are just initialized with a predefined value and these positions are not showing the actual locations. Laser scan is covering 210 degrees around the robot, so number of horizontal features is more than extracted vertical line features (which have only 120 degree view). The robot continues its motion. Figure 6(a) shows final location of the robot (having black square as its identity) while Fig. 6(b) illustrates 2D map with trajectory. The updated laser lines (in blue), laser lines corners (in cyan), vertical image lines showing inverse-depth points (in green) and current robot pose (in red) can be seen in the figure. New vertical lines (in green) are also observed which were not detected at the beginning. Each inverse-depth point, initially initialized with a default value, has now converged to its nearest world location. Successive observations lead to depth estimation and points shift to almost correct locations. Some lines are not converged correctly due to low parallax and problems in data association. Figure 7 presents the resultant partial 3D updated map of the environment and trajectory within it as shown by red marks.   In second experiment, another similar robot is also deployed. Assuming that this robot is moving from a known location with a known orientation in the same environment, the first image view from the second robot is shown in Fig. 8(a) with the corresponding laser scan illustrated in Fig. 8(b). In the image, the first robot can also be seen in the view of the second robot (red mark).  Second robot continues its motion and keeps building its map. Due to known location of the second robot w.r.t. the first one, the second robot poses and its map can be plotted in the same map model of the first robot (shown in Fig. 7). Figure 9 shows trajectories of the first and second robot in red and blue respectively in the resultant partial 3D global map of the environment. It is important to note that vertical and horizontal lines are increased as second robot has explored environment with different views as compared to first robot exploration. Table I shows the mean features position error when compared with their actual locations in the environment with the developed map model.

VI. CONCLUSIONS
A unified EKF SLAM framework is proposed and tested in an indoor environment for building partial 3D map by using more than one robot. The proposed method provides a compact map model as compared to other representations. It is easy to process and to handle during computation and sharing. Based on experiments, results are found nearly accurate for both kinds of sensors extracted features. For better performance of the proposed method, it is noticed that for image features, extraction and association of more features are required. The developed map model finds its potential in various applications requiring real-time survey of environment, 3D modeling, building map for robot navigation and path planning. Results can also beneficiate human assistance to accomplish various tasks.


HEURISTIC TECHNIQUES FOR HANDWRITTEN SIGNATURE CLASSIFICATION

New theoretical and experimental techniques for offline classification of handwritten signatures are introduced in this paper. The proposed algorithms are mainly based on boundary tracing technique for extracting characteristic features. Outer and inner boundaries of the signature image are treated separately. The upper and lower parts of the boundaries are extracted to form two sequences of points. Three algorithms for calculating feature vectors are applied based on y coordinate, distances between consecutive points and from polar coordinates system. Experiments on classification of the resulted vectors were carried out by means of Dynamic Time Warping algorithm using window and slope constraints. A brief comparison between the authors' work and other known signature techniques is also discussed in the paper.

INTRODUCTION
The handwritten signature is still very common way for authorizing various kinds of documents. From legal contracts to payment bills they play an important role and are used on everyday basis. Signatures are usually written on regular paper without any specialized equipment. Therefore, the only information available, which future verification may be based on, is the static image of the signature. This kind of authorization is obviously far from being perfect. It is not difficult for a skilled person to forge someone's signature. The shape of the signature can be duplicated when one has access to the original signatures and enough time to train signing.

INPUT DATA AND PREPROCESSING
In order to prepare data for classification algorithm, the images of signatures are first stored as Portable Network Graphics files (Fig 1). This particular format for graphical files provides lossless compression that retains all important features without introducing distortions, and results in relatively small footprint. Images can be obtained by means of scanning devices from original documents. The segmentation of signatures from acquired scans is not considered in this work, but can be easily implemented by applying certain constraints on the position of the signature inside the analyzed document. Another problem is noise and defects caused by poor quality of documents and the scanning process. In our experiments we used threshold technique to eliminate minor distortions and convert images from grayscale into black-andwhite binary map.

DATA REDUCTION
The line of a signature in an image may consist of a large amount of pixels. Depending on the resolution and thickness of ink trace it can even reach a few dozen thousands of points. Classification of such a complex object may pose a very difficult task. In fact, most of the points don't give additional information and can be safely ignored. There are many techniques for reducing their number whilst preserving the most important  [2], view-based approach [3], and contour based techniques [4,5].
During thinning process most of the points comprising a particular object are removed to achieve one-pixel-width skeleton (Fig 2). This approach has many applications and has been widely used in cursive script recognition systems.

Fig. 2 -Signature (a) with its thinned version (b)
Projection based techniques capture the distribution of ink in an object by projecting its pixels onto different axes and summing their number or intensity values. The following figure (Fig. 3b) shows a projection of a signature image calculated by computing total number of pixels in every column of the picture.

Fig. 4 -Signature (a) with its upper (b) and down (c) views
Contour tracing algorithm follows the boundaries comprising object image and collects coordinates of their consecutive points (Fig 5a). By boundary we mean the group of the object's pixels that share at least one edge with the background of the image. In many cases a signature is made of several disconnected segments like letters, dashes, points, etc. Those different parts can be treated separately or can be concatenated to form one continuous object.
Sometimes it is also useful to separate outer boundaries from internal and consider upper and down parts of the boundaries separately (Fig 5b, 5c, 5d).

Fig. 5 -Examples of full contour (a), internal contour (b), upper concatenated contour (c), down concatenated contour (d)
Those transformations may reduce complexity of classification task by considering different components individually. In order to further reduce the number of points a simple sampling may be used by selecting every M-th value from the acquired sequences to form feature vectors (value of M denotes a step in the sampling process). In this work we focused on boundary tracing for reduction of data. Experiments with other techniques were carried out and presented in [4,6].

FEATURE VECTOR COMPUTATION
For the classification process each of the signatures is represented by one or more vectors. There are various methods for constructing such vectors. These methods should preserve all the features necessary for distinguishing between different classes of signatures. The other goal is to improve the separation of the classes by ignoring disturbances created by roughness of the ink trace and minor artifacts, which mostly are conducive to inaccuracy of the signing individual. By the class of a signature we mean the group of signatures signed by a particular person.

y-coordinates
The first approach presented in this work is collecting y coordinates of subsequent pixels that form the boundary of analyzed signatures (Fig 6).
As a result a vector V is obtained as a signature representation used for classification process (1): n n y y y y V , ,..., ,

Consecutive points
Another technique is to describe subsequent points as vectors. These vectors are computed as the difference between positions of each consecutive pair of points (2,3). This process is illustrated by Fig. 7. [ ]

Polar coordinates
The third alternative is the algorithm used in this work is the calculation of vectors describing points as vectors from the origin of the polar coordinate system (Fig 8).
In addition, two modifications were applied to reduce amount of unintuitive alignments called singularities [7]. The first used a window which constrained possible paths in the matrix of The second used a slope constraint allowing warping path to follow only particular directions. The applied slope constraint [7,8] can be expressed by the following equation (Eq. 6):

RESULTS
In order to evaluate effectiveness of presented methods several experiments were carried out. The database of signatures was created by 20 different people, with each signature repeated three times, giving a total of 60 signatures.
For each person, each two of the signature versions were used as reference patterns to classify the third one. Therefore 60 3 20 = × tests were conducted in each variant of the experiment. In all cases boundary tracing algorithm was applied to reduce dimensionality of data. During boundary tracing a sampling step of 10 was used to eliminate disturbances created by roughness of the ink trace and to further reduce amount of redundant information.
In the first three approaches only external boundaries were used. Upper parts of boundaries were concatenated and treated separately from bottom parts as described in section 3. The classification process was based on distance measure computed by means of Dynamic Time Warping algorithm. The distance of the reference vector from the one being classified was calculated according to Eq. 7: When comparing subsequent vectors representing points in method 2 and 3 a measure given by Eq. 8 was applied:  table 1 show the percentage of properly classified signatures using each of the methods described earlier. The best classification rate was achieved by using only y coordinates of selected boundary points.
In addition to tracing external boundaries, contours of internal elements like loops were also examined. During experiments it was discovered that most signatures used to construct database were written carelessly. Many of the signatures written by one person varied in number and shape of internal loops. However, some consistency in the numbers of loops in signatures was found. To include this information in a distance measure we added the following modification (9): The following table (Tab. 2) shows results for different value of k using the first approach for building feature vectors. As can be seen from Tab. 2, changing the value of the weighting factor k has noticeably changed the classification rate. It is 93% for . When individuals are required to write their signatures more carefully, the applied algorithms for comparing shapes and loop position definitely result in even better classification rates.

COMPARISON WITH PREVIOUS WORKS AND OTHER APPROACHES
In our previous works [11,12] we experimented with view-based approach and several other modifications of DTW algorithm. Table 3 presents the results achieved with different versions of DTW and various window sizes [7,8].
The column titles have the following meanings: 1. p -percentage of maximum window's size, 2. DTW -basic algorithm with window constraint, 3. DDTW -dynamic time warping with window constraint, 4. DTWS -basic algorithm with window and slope constraints. The results show that applying boundary tracing techniques to signature image gives better classification rate than view-based approach. The advantage of the view-based method is that it is computationally less expensive.

CONCLUSIONS
The methods presented in this paper are based on the techniques applied to cursive word recognition. The signature images are examined as whole words without segmentation into distinct letters or strokes. The aim of this experimental approach is to enable identification based on handwritten signature that would compare general appearance of the signatures. Most of the research in the area of automatic signature recognition is focused on verification and resistance to forgery. Those systems are based on features and techniques allowing for rejection of forged signatures. Combining comparison based on general appearance with more restrictive verification methods may result in more flexible systems capable of both identification and verification for different levels of requirements on exactness of signature repetitions.
The new experiments and comparison with our recent works and other approaches show that the results obtained in [13] using the described approach had really been promising as presented there. The classification rate is still higher than the previous work results as shown in Tables 2 and 3. The authors' current research is being focused on incorporating other classification methods like Neural Networks or Hidden Markov Models. Toeplitz matrix minimal eigenvalues are still under studying to consider their use in feature points extraction. In addition, it is planned to combine and fuse the offline information collected from the signature image with the online data and information obtained from a camera or tablet devices in a hybrid system. These approaches will definitely increase the recognition rate as they had already done with other applications [9,10].


extfieldsofstudy:
Mathematics
Physics
Computer Science

text:
Tensor calculus with open-source software: the SageManifolds project

The SageManifolds project aims at extending the mathematics software system Sage towards differential geometry and tensor calculus. Like Sage, SageManifolds is free, open-source and is based on the Python programming language. We discuss here some details of the implementation, which relies on Sage's parent/element framework, and present a concrete example of use.

Introduction
Computer algebra for general relativity (GR) has a long history, which started almost as soon as computer algebra itself in the 1960s. The first GR program was GEOM, written by J.G. Fletcher in 1965 [1]. Its main capability was to compute the Riemann tensor of a given metric. In 1969, R.A. d'Inverno developed ALAM (for Atlas Lisp Algebraic Manipulator ) and used it to compute the Riemann and Ricci tensors of the Bondi metric. According to [2], the original calculations took Bondi and collaborators 6 months to finish, while the computation with ALAM took 4 minutes and yielded the discovery of 6 errors in the original paper. Since then, numerous packages have been developed: the reader is referred to [3] for a review of computer algebra systems for GR prior to 2002, and to [4] for a more recent review focused on tensor calculus. It is also worth to point out the extensive list of tensor calculus packages maintained by J. M. Martin-Garcia at [5].
All applications listed in the second category are free software. In the first category, xAct and Ricci are also free software, but they require a proprietary product, the source code of which is closed (Mathematica).

The SageManifolds project 4.1. Aim and scope
Sage is well developed in many areas of mathematics but very little exists for differential geometry and tensor calculus. One may mention differential forms defined on a fixed coordinate patch, implemented by J. Vankerschaver [20], and the 2-dimensional parametrized surfaces of the 3dimensional Euclidean space recently added by M. Malakhaltsev and J. Vankerschaver [21]. The aim of SageManifolds [22] is to introduce smooth manifolds and tensor fields in Sage, with the following requirements: (i) one should be able to introduce various coordinate charts on a manifold, with the relevant transition maps; (ii) tensor fields must be manipulated as such and not through their components with respect to a specific (possibly coordinate) vector frame.
Concretely, the project amounts to creating new Python classes, such as Manifold, Chart, TensorField or Metric, to implement them within the parent/element pattern and to code mathematical operations as class methods. For instance the class Manifold, devoted to real smooth manifolds, is a parent class, i.e. it inherits from Sage's class Parent. On the other hand, the class devoted to manifold points, ManifoldPoint, is an element class and therefore inherits from Sage's class Element. This is illustrated by the inheritance diagram of Fig. 1  diagram, each class at the base of some arrow is a subclass (also called derived class) of the class at the arrowhead. Note however that the actual type of a parent is a dynamically generated class taking into account the mathematical category to which it belongs. For instance, the actual type of a smooth manifold is not Manifold, but a subclass of it named Manifold with category, reflecting the fact that Manifold is declared in the category of Sets 2 . Note also that the class Manifold inherits from Sage's class UniqueRepresentation, which ensures that there is a unique manifold instance for a given dimension and given name.

Implementation of charts
Given a smooth manifold M of dimension n, a coordinate chart on some open subset U ⊂ M is implemented in SageManifolds via the class Chart, whose main data is a n-tuple of Sage symbolic variables (x 1 , . . . , x n ), each of them representing a coordinate. In general, more than one (regular) chart is required to cover the entire manifold. For instance, at least 2 charts are necessary for the n-dimensional sphere S n (n ≥ 1) and the torus T 2 and 3 charts for the real projective plane RP 2 (see Fig. 6 below). Accordingly, SageManifolds allows for an arbitrary number of charts. To fully specify the manifold, one shall also provide the transition maps (changes of coordinates) on overlapping chart domains (SageManifolds class CoordChange).

Implementation of scalar fields
A scalar field on manifold M is a smooth mapping  where U is some open subset of M. A scalar field has different coordinate representations F , F , etc. in different charts X,X, etc. defined on U : These representations are stored in some attribute of the class ScalarField, namely a Python dictionary 3 whose keys are the various charts defined on U : f. express = X : F,X :F , . . . .
Each representation F is an instance of the class FunctionChart, which resembles Sage native symbolic functions, but involves automatic simplifications in all arithmetic operations. Given an open subset U ⊂ M, the set C ∞ (U ) of scalar fields defined on U has naturally the structure of a commutative algebra over R: it is clearly a vector space over R and it is endowed with a commutative ring structure by pointwise multiplication: The algebra C ∞ (U ) is implemented in SageManifolds via the parent class ScalarFieldAlgebra, in the category CommutativeAlgebras.
The corresponding element class is of course ScalarField (cf. Fig. 2).

Modules and free modules
Given an open subset U ⊂ M, the set X (U ) of all smooth vector fields defined on U has naturally the structure of a module over the algebra C ∞ (U ). Let us recall that a module is similar to a vector space, except that it is based on a ring (here C ∞ (U )) instead of a field (usually R or C in physical applications). Of course, every vector space is a module, since every field is a ring. There is an important difference though: every vector space has a basis (as a consequence of the axiom of choice), while a module does not necessarily have any. When it possesses one, it is called a free module. Moreover, if the module's base ring is commutative, it can be shown that all bases have the same cardinality, which is called the rank of the module (for vector spaces, which are free modules, the word dimension is used instead of rank ).
If X (U ) is a free module (examples are provided in Sec. 4.5 below), a basis of it is nothing but a vector frame (e a ) 1≤a≤n on U (often called a tetrad in the context of 4-dimensional GR): The rank of X (U ) is thus n, i.e. the manifold's dimension 4 . At any point p ∈ U , Eq. (5) gives birth to an identity in the tangent vector space T p M: which means that the set (e a (p)) 1≤a≤n is a basis of T p M. Note that if U is covered by a chart (x a ) 1≤a≤n , then (∂/∂x a ) 1≤a≤n is a vector frame on U , usually called coordinate frame or natural basis. Note also that, being a vector space over R, the tangent space T p M represents another kind of free module which occurs naturally in the current context. It turns out that so far only free modules with a distinguished basis were implemented in Sage. This means that, given a free module M of rank n, all calculations refer to a single basis of M . This amounts to identifying M with R n , where R is the ring over which M is based. This is unfortunately not sufficient for dealing with smooth manifolds in a coordinate-independent way. For instance, there is no canonical isomorphism between T p M and R n when no coordinate system is privileged in the neighborhood of p. Therefore we have started a pure algebraic part of SageManifolds to implement generic free modules, with an arbitrary number of bases, none of them being distinguished. This resulted in (i) the parent class FiniteRankFreeModule, within Sage's category Modules, and (ii) the element class FiniteRankFreeModuleElement. Then both classes VectorFieldFreeModule (for X (U ), when it is a free module) and TangentSpace (for T p M) inherit from FiniteRankFreeModule (see Fig. 3).

Implementation of vector fields
Ultimately, in SageManifolds, vector fields are described by their components with respect to various vector frames, according to Eq. (5), but without any vector frame being privileged, leaving the freedom to select one to the user, as well as to change coordinates. A key point is that not every manifold admits a global vector frame. A manifold M, or more generally an open subset U ⊂ M, that admits a global vector frame is called parallelizable. Equivalently, M is parallelizable if, and only if, X (M) is a free module. In terms of tangent bundles, parallelizable manifolds are those for which the tangent bundle is trivial: Examples of parallelizable manifolds are [23] • the Cartesian space R n for n = 1, 2, . . ., • the sphere S 3 SU(2), as any Lie group, • the sphere S 7 , • any orientable 3-manifold (Steenrod theorem [24]).  • the sphere S 2 (as a consequence of the hairy ball theorem), as well as any sphere S n with n ∈ {1, 3, 7}, • the real projective plane RP 2 .
Actually, "most" manifolds are non-parallelizable. As noticed above, if a manifold is covered by a single chart, it is parallelizable (the prototype being R n ). But the reverse is not true: S 1 and T 2 are parallelizable and require at least two charts to cover them.
If the manifold M is not parallelizable, we assume that it can be covered by a finite number In particular, this holds if M is compact, for any compact manifold admits a finite atlas. We then consider the restrictions of vector fields to the U i 's. For each i, X (U i ) is a free module of rank n = dim M and is implemented in SageManifolds as an instance of VectorFieldFreeModule (cf. Sec 4.6. Implementation of tensor fields The implementation of tensor fields in SageManifolds follows the strategy adopted for vector fields. Consider for instance a tensor field T of type (1,1) on the manifold M. It can which defines T a b , is meaningful only when a vector frame (e a ) exists 5 . Therefore, one first decomposes the tensor field T into its restrictions T | U i on parallelizable open subsets of M, U i (1 ≤ i ≤ N ) and then considers the components on various vector frames on each subset U i . For each vector frame (e a ), the set of components (T a b ) is stored in a devoted class (named Components), which takes into account all the tensor monoterm symmetries: only non-redundant components are stored, the other ones being deduced by (anti)symmetry. This is illustrated in Fig. 5, which depicts the internal storage of tensor fields in SageManifolds. Note that each component T a b is a scalar field on U i , according to the formula Accordingly, the penultimate level of Fig. 5 corresponds to the scalar field storage, as described by (3). The last level is constituted by Sage's symbolic expressions (class Expression).

Parallelization
To improve the reactivity of SageManifolds and take advantage of multicore processors, some tensorial operations are performed by parallel processes. The parallelization is implemented by means of the Python library multiprocessing, via the built-in Sage decorator @parallel. Using it permits to define a function that is run on different sub-processes. If n processes are used, given a function and a list of arguments for it, any process will call the function with an element of the list, one at time, spanning all the list.
Currently 6 , the parallelized operations are tensor algebra, tensor contractions, computation of the connection coefficient and computation of Riemann tensor.
The parallelization of an operation is achieved by first creating a function which computes the required operation on a subset of the components of a tensor; second, by creating a list of 2n (twice the number of used processes) arguments for this function. Then applying this function to the input list, the calculation is performed in parallel. At the end of the computation a fourth phase is needed to retrieve the results. The choice to divide the work in 2n is a compromise between the load balancing and the cost of creating multiple processes. The number of processors to be used in the parallelization can be controlled by the user.

SageManifolds at work: Kerr spacetime and Simon-Mars tensor
We give hereafter a short illustration of SageManifolds focused on tensor calculus in 4-dimensional GR. Another example, to be found at [25], is based on the manifold S 2 and focuses more on the use of multiple charts and on the treatment of non-parallelizable manifolds. Yet another example illustrates some graphical capabilities of SageManifolds: Figure 6 shows the famous immersion of the real projective plane RP 2 into the Euclidean space R 3 known as the Boy surface. This figure has been obtained by means of the method plot() applied to three coordinate charts covering RP 2 , the definition of which is related to the interpretation of RP 2 as the set of straight lines ∆ through the origin of R 3 : (i) in red, the chart X 1 covering the open subset of RP 2 defined by all lines ∆ that are not parallel to the plane z = 0, the coordinates of X 1 being the coordinates (x, y) of the intersection of the considered line ∆ with the plane z = 1; (ii) in green, the chart X 2 covering the open subset defined by all lines ∆ that are not parallel to the plane x = 0, the coordinates of X 2 being the coordinates (y, z) of intersection with the plane x = 1; (iii) in blue, the chart X 3 covering the open subset defined by all lines ∆ that are not parallel to the plane y = 0, the coordinates of X 2 being the coordinates (z, x) of intersection with the plane y = 1. Figure 6 actually shows the coordinate grids of these three charts through the Apéry map [26], which realizes an immersion of RP 2 into R 3 . This example, as many others, can be found at [25].

function(arguments)
In a functional language, this would correspond to result = function(object,arguments). For instance, the Riemann tensor of a metric g is obtained as riem = g.riemann() (in this case, there is no extra argument, hence the empty parentheses). With this in mind, let us proceed with the computation by means of SageManifolds. In the text below, the blue color denotes the outputs as they appear in the Sage notebook (note that all outputs are automatically L A T E Xformatted by Sage).
The first step is to declare the Kerr spacetime (or more precisely the part of the Kerr spacetime covered by Boyer-Lindquist coordinates) as a 4-dimensional manifold:
Let us conclude by stating that, in the very spirit of free software, anybody interested in contributing to the project is very welcome!


extfieldsofstudy:
Computer Science
Medicine

text:
Feature Optimization for Long-Range Visual Homing in Changing Environments

This paper introduces a feature optimization method for robot long-range feature-based visual homing in changing environments. To cope with the changing environmental appearance, the optimization procedure is introduced to distinguish the most relevant features for feature-based visual homing, including the spatial distribution, selection and updating. In the previous research on feature-based visual homing, less effort has been spent on the way to improve the feature distribution to get uniformly distributed features, which are closely related to homing performance. This paper presents a modified feature extraction algorithm to decrease the influence of anisotropic feature distribution. In addition, the feature selection and updating mechanisms, which have hardly drawn any attention in the domain of feature-based visual homing, are crucial in improving homing accuracy and in maintaining the representation of changing environments. To verify the feasibility of the proposal, several comprehensive evaluations are conducted. The results indicate that the feature optimization method can find optimal feature sets for feature-based visual homing, and adapt the appearance representation to the changing environments as well.

Introduction
A usual dichotomy partitions robot visual navigation methods into quantitative and qualitative methods [1]. The estimation of metric position information is not needed in a qualitative method, and the environments are usually represented by a route. It is crucial for robots to move autonomously from their current position to the home position along the route. Numerous methods for route guidance have been proposed, including node-based [2][3][4][5][6] ones and those that do not use nodes [7,8]. Navigation models that do not separate routes into sequences of waypoints have been proposed in [7,8], following the widely neglected branch of models based on recognition-triggered response. The navigation with nodes has been modeled by local visual homing strategies [2][3][4][5]. Inspired by the simple principle of biological navigation, local visual homing has the ability to return to a previously visited location, and has drawn attention in robotics study as computationally cheap building blocks in a minimalistic appearance-based navigation framework.
In this paper, we present that long-range homing along a route can be achieved by implementing local visual homing strategy to sequentially visit the topological nodes. In the feature-based local visual homing in unprepared environments, the extracted visual features are ideal alternatives to artificial landmarks. However, a big challenge during long-range homing in changing environments is the management of rich feature sources [9,10]. In order to avoid irrelevant and redundant features, task-oriented selection schemes are needed, which contain additional constraints in evaluating feature usefulness for the tasks. Furthermore, in order to obtain optimal features during long-range navigation, it is necessary to maintain an up to date representation of surrounding environments. Propelled by above requirements, the main work of this paper focuses on the feature optimization orienting towards long-range feature-based visual homing in changing environments.
Several feature-based visual homing strategies have been presented in the literature, such as average displacement vector [11], average landmark vector (ALV) [12] and the included angle difference method [2]. However, they all need or imply an equal-distance assumption that all landmarks are distributed at the same distance from the snapshot location. Since the assumption is always violated in unstructured environments, the true home vector will deviate from the estimated direction [2,11,13]. Uniformly distributed features do improve the feature-based homing accuracy, especially for the ALV, which has drawn significant attention for simplicity [13]. There are adequate visual features extracted from unstructured environments to identify a goal to home, yet no attention has been given to modify the feature distribution for feature-based visual homing. A few of the widely employed feature extraction algorithms have tried to get uniformly distributed features [14][15][16], however, none meets the requirement on feature distribution for feature-based visual homing.
The feature selection is essential to object classification [17], localization [18], robot navigation [19,20], etc. With respect to the features used in visual homing, there are several criteria for the selection process [21,22], in which the demanding task is the explicit quantitative characterization of feature properties in view of their relative importance. In the previous literature, most approaches to the characterization of feature quality focus on recognition and classification tasks, but few of them are ideally suited to feature-based visual homing. Furthermore, most previous research on visual homing was carried out under the assumption that the environments are static. The significance of rating and updating mechanisms, which are crucial to continuously evaluating the relative importance of features and discarding useless ones, is always ignored.
Motivated by the aforementioned thought, our work concerns the optimization of feature distribution, selection and updating. In particular, we focus on acquiring uniformly distributed features to fulfill the equal-distance assumption. The features are graded by the quantitatively characterized selection criteria of visual homing. When the agent retraces the environments, the importance of features is re-evaluated to update the appearance representation. In this paper, the ALV strategy is adopted as building blocks of the route for simplicity. Besides, the features are extracted by SURF algorithm [23], because of its high accuracy and less computing time. In order to improve the performance of long-range feature-based visual homing in changing environments, the work presented in this paper concentrates on maximizing the advantage of the ALV method by modifying the distribution of high quality SURF features.
The remainder of the paper is organized as follows: Section 2 outlines the extraction algorithm of well-distributed SURF features. Section 3 presents the feature selection and updating mechanisms. Section 4 shows the framework of feature optimization. Experiments in Section 5 demonstrate the performance. Section 6 draws conclusions and points out future work directions.

Uniformly Distributed Features
Due to their good invariance, local features have been introduced to substitute for artificial landmarks when the agent is situated in unknown environments. The Speeded Up Robust Features (SURF) [23] and Scale-Invariant Feature Transform (SIFT) [24] are widely used, however, neither controls the feature distribution. Several modified versions have focused on improving spatial distribution [14][15][16], but none can ideally meet the need of feature-based visual homing. In addition, little work has been done to improve the distribution of SURF features. Considering that SURF burden is computationally low, we aim at improving the spatial distribution. Considering that the SURF features are extracted in scale space and the image sizes are identical, extracting features over the scale space in a uniform way becomes very important. The SURF scale space consists of several octaves and scale layers, as shown in the left side of Figure 1.
In consideration of using panoramic images, each scale layer can be divided into regular sector rings along the radial and circumferential directions (right side of Figure 1). The features should be uniformly distributed in each sector ring. Therefore, the feature distribution problem in the image can be seen as a problem in scale layers. Since the scale space is constructed by up-scaling the box filter size, the number of features decreases as the filter size increases because of the smoothing characteristics. Besides, the number of sector rings decreases progressively among octaves of the scale space.
Supposing that the number of key-points extracted by the standard SURF algorithm in the scale layer (s) of the octave (o) is N os , the number of sector rings is n os , and the key-point number in the ith sector ring is M i os . The feature distribution in the scale layer will be relatively uniform if the key-point number in each sector ring is the same. We can define N i os in Equation (1), which denotes the target number of key-points when each sector ring has the same key-point number. For the ith sector ring, all the key-points in this sector ring will be reserved when M i os is not greater than N i os . On the contrary, if M i os is greater than N i os , the excess key-points will be discarded by their quality which is measured according to strength value and spatial dispersion. The strength value V str will be shown in Equation (16) and specifically explained in Section 3.1. With regard to spatial dispersion, it is computed as follows. As shown in Equation (2), E j is the entropy of the square region used to construct the descriptor of the jth key-point, and q l is the probability of the lth gray level value over all the grey levels within the square region: Supposing (x j , y j ) is the position of the jth key-point in the ith sector ring. The center ( , ) xy of all the key-points in the ith sector ring can be computed by Equation (3) (3) where E j is taken as the weight of the jth key-point. Supposing dist j in Equation (4) is the Euclidean distance between the jth key-point and the center, and Dist in Equation (5) is the divergence factor of the ith sector ring. Then the quality of the jth key-point in the ith sector ring can be measured by its strength value and spatial dispersion, as given in Equation (6) where W D is weight factor: In order to acquire uniformly distributed key-points, the points with lower quality are rejected to make sure that the key-point number is less than N i os in each sector ring. A pseudo-code description of the extraction of uniformly distributed features is given in Algorithm 1. if TRUE then 8: Calculate the key-point quality according to Equation (6). 9: Classify the key-points by their quality in descending order.

Proposed Approaches to Feature Selection and Updating
The feature strength is crucial in feature-based visual homing in unstructured environments where adequate visual features can be extracted. With respect to the features used in the homing process, there should be several criteria of feature selection to evaluate their strength, in addition to the feature distribution. Besides, the fulfillment of long-range navigation in changing environments critically depends on the feature updating.

Feature Selection
Besides the criterion mentioned above, there are three more criteria identified for selection of features [21,22]: (1) relevance: the relevance of a feature is defined as its importance in making navigational decisions; (2) distinctiveness: features should be clearly distinguishable from their surroundings, which means the feature configurations should be unique; (3) reliability: taking the invariance of feature descriptor into account, features should cope with normal image deformations, which requires high robustness. What really counts is that the features should be widely visible, and consistently observed when the robot passes the same position every time.
The next step is to get the quantitative characterization of the criteria, in order to assign individual strength value to each feature. From the relevance aspect, individual feature evaluation is based on their relevance to the task and relationships with other features. Relevance and relationships are generally characterized in terms of correlation and mutual information. The distribution quality, discussed in Section 2, is significant for visual homing and has been taken into account during feature selection process, thus every selected feature is surely relevant to the homing task. The mutual information has been applied for navigation tasks [6,20,25,26], however, it is used in a different way in this paper. The mutual information can be used to measure the quantity of information shared by two random variables. Therefore, it can be used to evaluate the relationships among features. The more information it measures, the more relevant the feature relationships are. The mutual information is defined as follows: where the entropy H(X) measures the variability of variable X. H(XY)is the conditional entropy, giving the expected information in X if a measurement is taken.
Supposing C denotes the group of all the T features in one sector ring. The mutual information MI(C,F) between group C and feature F in this group can be measured by Equation (7). Due to the lack of previous knowledge, it is supposed that p(C) has a uniform prior. Therefore, all the T features in group C have the same probability and the entropy of group C can be computed by Equation (8). In addition, the conditional entropy is given by Equation (9): where p(CF) can be computed by the Bayes rule: Here, p(CF) is estimated by applying a kernel density estimator [27]: In Equation (11), K(•) is the Gaussian Kernel function with standard deviation, and d(•) denotes the Euclidean distance. An efficient and accurate nearest-neighbor approximation of this estimator is done as follows [27]: During the navigation, the images captured at current positions are continuously matched with the image stored at the target. Supposing M is the number of times the target image has been matched, and m is the number of times that the feature F extracted from the target image has been matched. The mutual information of feature F can be computed according to Equation (13): where F j is the matching feature of F, d j denotes the matching distance between F and F j , and the structural principle of C j is the same as C.
A distinctiveness selection can be on the basis that the more informative the feature is, the greater its distinctive power reduces. The feature distinctiveness is written as Equation (14): The contrast value indicates the reliability property directly. In addition, the power of reliability can also be estimated by the record times that the feature has been successfully re-observed. Thus the quantitative characterization of reliability can be measured according to Equation (15): where H norm denotes the normalized Hessian value (SURF key-point), and W H is the weight factor. The feature strength value V str , given in Equation (16), is obtained by adding the three quantitative characterizations: where MI norm is the normalized mutual information, D norm is the normalized reliability, W M and W R are the mutual information and the reliability weight factors, respectively.

Feature Updating
As for the long-range navigation in changing environments, the feature updating is a matter of cardinal significance. The approaches close to our work are presented in [9,10]. They are both based on the multi-store model of human memory [28], which is divided into three basic stores: sensory memory (SM), short-term memory (STM) and long-term memory (LTM). The sensory memory stores all information perceived by the senses, and elementary selection measures are taken to determine what information will be transported to the STM. Through a process of rehearsal, the information in STM can be transported to LTM to be retained for longer periods of time. The approach presented in [9] is mainly used to keep feature sets updating, and the feature upgrade or degradation is based on feature "hit" or "miss" rather than quantitative characterization. An adapted memory model is adopted to update the features in [10], where the main purpose is robot localization.
In this paper, the feature updating is implemented by adopting the concepts mentioned above in a simple way. During the process of SURF feature extraction, original key-points are denoted as SM, and those with low contrast and large principal curvature will be rejected. Then the final features are classified as STM or LTM according to their strength values. For the purpose of finding the most suitable threshold T str to distinguish STM from LTM, an experiment described specifically in Section 5.1 is conducted to test the homing performance. When the feature strength values are re-computed, the feature sets will be updated. For STM, they have four times to be transported to LTM, on the premise that their strength values are greater than the threshold T str , otherwise, they will be deleted. Besides, the LTM will degrade to STM if the strength values are less than the threshold T str .

Overview of the Proposed Feature Optimization Procedure
The outline of the feature optimization procedure of long-range visual homing is shown in Figure 2. The route consisting of nodes is determined when the robot enters the environments for the first time. When the robot retraces the environments, it follows the route using the strategies in [2][3][4][5][6]. During the retracement process, the ALV homing algorithm directs the robot towards the target node. In addition, the ALV homing direction is computed with the assistance from an external compass. The two-stage procedure is detailed as follows: (1) Map building. The robot is manually driven in the environments and panoramic images are captured at particular places with approximately 2 m between each other. Each node in the route is a collection of optimal features, which are extracted from the panoramic image. In the first step, each panoramic image is divided into regular sector rings as shown in Figure 1, and features are extracted by standard SURF algorithm. Then the strength value is calculated for each feature, and all features are classified into STM and LTM according to their strength values.
For all the LTM features, the method proposed in Section 2 is used to improve the feature distribution, regardless of the STM features. If the LTM feature number in one sector ring is less than the number required, the STM features with higher strength values compensate for the shortage of feature number. The uniformly distributed LTM features are used to compute the ALV vectors.
(2) Retrace phase. When the robot follows the succession of nodes by ALV homing method, the feature updating approach will be adopted.
During the navigation between two nodes, panoramic images are continuously captured and matched with the image taken at the target node before the robot arrives at the node. The strength values of features stored at the target node are re-evaluated during this process. Therefore, the feature set can be updated. It should be noticed that feature matching aims at re-evaluating the strength values rather than calculating homing direction, since there is no need of feature matching in the ALV method.
The homing direction will be changed in two cases. One case is that, if the LTM features used to compute current direction of migration change by T num or higher, the ALV vector at current location will be computed to obtain new homing direction. The other case is that, if the LTM features stored at the target node is updated, its ALV vector will be re-computed and the robot will move on in a new direction.
With regard to route guidance using linked local visual homing, there are two issues of crucial importance that have to be noticed: the determination of a new node and the arrival detection at a node [3,7,8]. As part of autonomous navigation, the major objective of this work is to maintain feature optimization, thus, the determinations on at which point a new node is set and whether the robot is sufficiently close to the target node are not the key points. In this paper, the locations of nodes are manually selected at a distance of approximately 2 m. Besides, the arrival detection of nodes is straight. During navigation between nodes, features extracted at current position are continuously matched with features stored at the target node. Arrival is declared once the maximum value of the bearing angle differences of matching features falls below threshold T θ .

Robot Platform and Parameter Settings for Experiments
Our experiments were performed with a crawler mobile robot (see Figure 3). A panoramic vision system, composed of a mirror with a diameter of 100 mm and a camera with a resolution of 800 × 800, was mounted on top of the robot. A 2 GHz Pentium (R) processor was dedicated to vision processing and robot control. In addition, the robot was equipped with a digital compass to acquire the global orientation.  Table 1 shows the parameter settings. Similar to the standard SURF algorithm, the scale space consists of three octaves and four scale layers (O = 3 and S = 4). Because the number of features decreases gradually with the number of octaves, the number of sector rings n os should taper off. More specifically, the sector ring number in each scale layer of the first octave is 100 (5 × 5 × 4), and the others are 64 (4 × 4 × 4) and 36 (3 × 3 × 4). The weight factor W D in Equation (6) must be determined by an appropriate compromise between the spatial dispersion and feature strength value. The optimal value of W D was determined by a local homing experiment, where we used 60 random pairs of panoramic images captured in the room and corridor. When the homing error was calculated, the Hessian value was substituted for feature strength value in Equation (6). The value assigned to W D was changed from 0 to 1 in 0.1 step. By seeking the minimum average angular error of homing, 0.3 was determined as the optimal value. The same experiment was conducted to determine the optimal value of W H , and the only difference was that the reliability value was substituted for feature strength value in Equation (6). Finally W H was set to 0.6 through the experimental analysis. Similar to the above two experiments, another experiment was performed to figure out the most suitable W M and W R values. The feature strength value measured according to Equation (16) was used. The minimum average angular error of homing was obtained when W M value was 0.5 and W R value was 0.3. We conducted another experiment to find the optimal threshold to distinguish LTM features from STM features. The threshold T str must balance feature quantity and quality. A low threshold will result in redundant features, and the features with low relative importance will degrade the homing performance. On the contrary, a high threshold will prune too many features, giving rise to the lack of representativeness of environmental appearance. The STM features will be less accessible to turn into LTM ones. The experimental results are shown in Figure 4. A well-established compromise in the mutual contradiction mentioned earlier was reached when the threshold T str was set to 0.7. In addition, the thresholds T num for changing homing direction and T θ for arrival detection were empirically set to 10% and 10°, respectively.

Evaluation Tests of the Uniformly Distributed Features
To evaluate the improvement on feature distribution using the proposed approach, an experiment was conducted in the room with a 10 m × 5 m floor space and the corridor with a 16 m × 3 m floor space. Figure 5a,c show the distribution of features extracted by the standard SURF algorithm. They manifest that the feature distribution is non-uniform, which means that the homing performance will be degraded. Figure 5b,d show the feature extraction results by the proposed modified version of SURF algorithm. The comparison shows the feasibility of the proposed approach in the improvement on feature distribution. However, this is just visual inspection. The improvement on the distribution quality should be judged in terms of visual homing. The chosen measure of quantitative characterization of distribution, similar to the one proposed in [29], is defined as:  (17) where p i is the ratio of matching features in the ith sector ring over all the ones in the scale layer (s) of the octave (o). If the correctly matching features were uniformly distributed, a maximum distribution quality would be achieved (5.94 for our experiments). (c) (d) Figure 6 shows the homing error and the corresponding feature distribution quality. Homing error is the absolute angular error between correct homing vector and the computed homing vector. In addition, the homing method presented in [4] was taken for comparison, as shown in Figure 6e,f. The proposed approach performed as well as the method in [4], with slight performance differences. The promising results indicate that it is feasible to improve the homing accuracy of ALV by employing the proposed approach, which can decrease the impact of anisotropic feature distribution.  . (a,b) show the homing error in the room using features extracted by the standard SURF and the proposed approach, respectively; (c,d) show the distribution quality corresponding to the results in (a,b); (e,f) show the homing error in the room and the corridor, respectively using the method in [4]; (g,h) show the homing error in the corridor using features extracted by the standard SURF and the proposed approach, respectively; (i,j) show the distribution quality corresponding to (g,h).

Long-Range Homing Experiments
In this section, long-range homing experiments were conducted to demonstrate the feasibility of the feature optimization approach, and the standard SURF algorithm was taken for comparison. Figure 7 shows the environment representation marked with four typical regions, where no special arrangements and modifications were made to facilitate the execution of visual homing. The robot was located upon the starting node, and sequentially passed the intermediate nodes until it arrived at the home node eventually. Notice that the initial orientations of the robot were roughly the same in the following experiments. Figure 7. The experimental environments. Phase 1 denotes the visual homing process in the room and the junction between the room and corridor; phase 2 corresponds to the process in the junction between two corridors; phase 3 denotes the process in the corridor; phase 4 corresponds to the process in the hall.
The experiments were classified into several groups according to the degree of change of the environmental appearance. The environments were manually changed in various manners, including modifying illumination conditions by control of curtains and lamps, rearranging object locations, adding new objects, removing existing objects and covering objects with cardboard. It should be noticed that temporary occlusions caused by people always existed. Here, five types of experiments had been conducted: in the first one, under good illumination conditions, the environmental layout remained the same as that in the first tour without modification; the second one was conducted under worse illumination conditions, compared with the first one; in the third and fourth experiments, environment changes were identical, and the illumination conditions were the same as that in the first and second, respectively; the fifth one was conducted in randomly changing environments. By means of the proposed feature optimization approach, the features stored in nodes updated five times via the five types of experiments mentioned above.
In order to demonstrate the effectiveness of each part of the feature optimization procedure, our approach was compared with three other approaches. The experimental results are presented in Figures 8-12. The results of the proposed approach and the standard SURF are shown by red and violet, respectively. The blue indicates the results obtained by our approach without the procedure of feature selection and updating. And the results shown by green are obtained by our approach without the procedure of feature updating. The decrease of curvature and the improvement on smoothness can be chosen as criteria to the effectiveness of the proposed approach on homing accuracy. Figures 8 and 9 show the homing trajectories under different illumination conditions in non-modified environments. Figure 8 corresponds to the results obtained at around noon with fairly good illumination, and Figure 9 corresponds to the results obtained in the afternoon with relatively poor illumination. It can be observed that the experimental results obtained under fairly good illumination condition are better than those obtained under relatively poor illumination conditions.

Long-Range Homing with Different Illumination Conditions in Non-Modified Environments
Observing the results of the four approaches, our approach outperforms the other three ones. The performance of the standard SURF is the worst, which may be caused by its poor control of the feature distribution and quality. The performance differences between the two experimental results are relatively bigger when using the standard SURF, compared with the other three approaches. This can be attributed to the utilization of feature selection and updating, which partially decrease the influence of large changes in illumination and people. Moreover, it is noteworthy that the homing performance in the corridor is improved apparently owing to the uniformly distributed features, compared to that in the room and hall. Although the performance is improved in the corridor, the homing error is bigger compared with that in the room and hall. Besides, a problem has been arisen that there exists significant degradation of homing performance in the junctions of two corridors, between the room and corridor, since the environmental appearance information between them is obviously different. (c) (d)

Long-Range Homing with Different Illumination Conditions in Modified Environments
The experimental results under different illumination conditions in modified environments are provided in Figures 10 and 11. Figure 10 indicates the results obtained at around noon with fairly good illumination, and Figure 11 is obtained at around afternoon with relatively poor illumination. The modifications made in the environments will result in decreasing of matching feature number, thus the ability of coping with the large changes in environments can be tested. Analyzing the homing trajectories, our approach performs well despite the changes in environmental layout. Though the homing error caused by anisotropic feature distribution decreases when the robot approaches the target node, the performance differences demonstrate that the feature optimization procedure does take effect. In a word, our approach maintains a better appearance representation of the changing environments.  Figure 12 shows the homing trajectories under randomly changing environmental conditions. The illumination conditions were changed randomly by controlling the curtains and lamps, and the object positions were in disorder, posing a challenge to visual homing. Comparing the experimental results with those in Sections 5.3.1 and 5.3.2, it is clear that the performance differences increase as the feature set updating times goes up. It is worth to note that the approach using standard SURF is incapable of dealing with these situations. From the experimental results, the effect of the feature optimization in our approach, which is the result of these performance differences, is evident. Given these challenging conditions, our approach performs well and demonstrates that it is of importance to get a suitable appearance representation of the environments.

Analysis of Experimental Results
From the perspective of route following, the homing trajectories in this paper are not the straightest, compared with those in [30][31][32]. Considering the simple navigation framework, the homing error tends to decrease as the robot is guided closer to the target node. Besides, the method of selecting node locations may result in the lack of features in common between consecutive nodes. Therefore, the homing direction will be changed frequently due to the strategy of changing mobile direction. The changing environments, the extreme unevenness of feature distribution in the corridor, and the transition from one environment to another with obviously different appearance pose severe challenges to visual homing.
To quantify the performance, the mean values and standard deviations of the mobile direction angles during the long-range homing are listed in Table 2, where the smoother trajectory and better performance is indicated by lower mean value and standard deviation. From the homing trajectories and mobile direction angles, the feature optimization method improves the homing accuracy. However, it is hard to know how well the proposed method performs due to the error introduced by the digital compass. During the first tour when the route was determined, the images were captured when the robot was parallel to one side of the oblong floor tile. Therefore, the directions of translation between images stored at neighboring nodes are known. The homing errors between neighboring nodes can be computed using the feature sets obtained through the above five experiments, without impact of the digital compass. The results in Table 3 indicate that each part of the optimization procedure has a positive effect on the decrease of homing error. Table 2. Mobile direction angles. Experiments 1-5 correspond to the results shown in Figures 8-12, respectively. E (°) denotes the mean value of mobile direction angle, and σ (°) denotes the standard deviation of mobile direction angle. Compared with the standard SURF (violet trajectory), P E and P σ denote the percentages by which mean value and standard deviation have been lowered through the five experiments, respectively. Despite the homing performance improvement, the equal-distance assumption is difficult to meet completely. This is obvious in the corridor and can be explained in the way that features are projected on the panoramic images [13]. In addition, the mechanisms of arrival detection and new node determination in this paper do not perform very well. The robot does wrongly declare the node arrival in a number of cases. The set distance of 2 m between consecutive nodes does not work in every environment and results in lack of autonomy. Table 3. Average homing errors. Feature sets 1-5 correspond to the features stored at the nodes and obtained during the above five types of experiments. Compared with the standard SURF (violet trajectory), P denotes the percentage by which the homing error has been lowered through the five feature sets.

Conclusions
This paper proposes an approach to feature optimization, especially for long-range feature-based homing in changing environments. The task-dependent optimization procedure consists of feature distribution, selection and updating. The method proposed in this paper shows good homing performance, including improvement on feature distribution quality, homing accuracy and environmental adaptability. Therefore, it is feasible and suitable for long-range navigation in changing environments.
The selection of key images is significant for representing a route with fewer nodes. In addition to that, the arrival detection is closely related to homing accuracy. However, the locations of nodes are determined manually and arrival detection mechanism does not work well in this paper. As crucial parts of autonomous navigation, more efforts need to be spent on solving these two issues. Other possible improvements lie in computing the most suitable strength value threshold with an adaptive algorithm to distinguish STM features from LTM features, as well as exploring the possibility of substituting visual compass for inaccurate digital compass.