Localization of Mobile Robots Using Odometry and an External Vision Sensor

This paper presents a sensor system for robot localization based on the information obtained from a single camera attached in a fixed place external to the robot. Our approach firstly obtains the 3D geometrical model of the robot based on the projection of its natural appearance in the camera while the robot performs an initialization trajectory. This paper proposes a structure-from-motion solution that uses the odometry sensors inside the robot as a metric reference. Secondly, an online localization method based on a sequential Bayesian inference is proposed, which uses the geometrical model of the robot as a link between image measurements and pose estimation. The online approach is resistant to hard occlusions and the experimental setup proposed in this paper shows its effectiveness in real situations. The proposed approach has many applications in both the industrial and service robot fields.


Introduction
This paper presents a sensor system composed of a single camera attached to a fixed position and orientation in a bounded environment ( indoor workplace) which is observing a mobile robot. The aim of the sensor system is to obtain the orientation and position (i.e., pose) of a mobile robot using both visual information retrieved by the camera and relative odometry readings obtained from the internal sensors of the robot. The camera acquisition and image processing tasks are executed in a specialized hardware, which also controls the behavior and internal sensors of the mobile robot through a wireless channel. The proposed schema allows the robot to perform complex tasks without requiring dedicated processing hardware on it. This approach is sustained in the Intelligent Space [1] concept and it can be equally applied to multiple scenarios, specially both in the industrial field (e.g., automatic crane positioning, autonomous car parking) and in the service fields (e.g., wheelchair positioning in medical environments or autonomous driving of mobile platforms ). The single camera solution presented in this paper allows to cover large areas with less cameras compared to multiple camera configurations where overlapped areas are mandatory. This feature reduces the cost and improves the reliability of the intelligent space philosophy.
In this paper, we suppose that the camera is correctly calibrated and thus the parameters governing the projection model of the camera are previously known. To connect the pose of the robot with information found in the image plane of the camera, we propose to obtain a 3D geometric model of the mobile robot. Such model is composed of several sparse points whose coordinates represent some well-localized points belonging to the physical structure of the robot. These points are determined by image measurements, called image features, which usually correspond to corner-like points due to texture changes or geometry changes such as 3D vertexes. Usually in industrial fields the image features are obtained by including some kind of artificial marker on the structure of the robot ( infrared markers or color bands). These methods are very robust and can be used to recognize human activity and models with high degrees of freedom (AICON 3D online [2] or ViconPeak online systems [3]). However, in this paper we want to minimize the required "a priori" knowledge of the robot, so that it is not necessary to place artificial markers or beacons on it to detect its structure in the images.
Independent of the nature of the image features, the information obtained from a camera is naturally ambiguous and thus some sort of extra metric information are required in order to solve such ambiguity. In this work, we propose to use the odometry sensors inside the robot (i.e., wheel velocities) to act as the metric reference.
The general diagram of the algorithm proposed in this paper is shown in Figure 1. It shows a clear division of the processes involved in obtaining the pose of the robot: first we denote as "Initialization of Pose and Geometry" to those processes necessary to start up the system, such as the 3D model of the robot and the initial pose it occupies. The initialization consists of a batch processing algorithm where the robot is commanded to follow a certain trajectory so that the camera is able to track some points of the robot's structure under different viewpoints jointly with the recording of the odometry information. All this information is combined to give the 3D model of the robot and the initial pose it occupies.
Given the initialization information, the second group of processes, named "Sequential Localization", provides the pose of the robot in a sequential manner. It is composed of a pose estimator, given odometry readings and a pose correction block which combines the estimation of the pose with image measurements to accurately give a coherent pose with the measurements. This algorithm operates entirely on-line and thus the pose is available at each time sample.
Both group of processes are supplied with two main sources of information: 1. Image measurements: they consist of the projection in the camera's image plane of certain points of the robot's 3D structure. The measurement process is in charge of searching coherent correspondences through images with different perspective changes due to the movement of the robot.
2. Motion estimation of the robot: The odometry sensors built on-board the robot supply the localization system with an accurate motion estimation in short trajectories but that is prone to accumulative errors in large ones. Figure 1. General diagram of the proposed localization system using a vision sensor and odometry readings.

Previous Works
Despite the inherent potential of using external cameras to localize robots, there are relative few attempts to solve it compared to the approach that consider the camera on-board the robot [4,5]. However, some examples of robot localization with cameras can be found in the literature, where the robot is equipped with artificial landmarks, either active [6,7] or passive ones [8,9]. In other works a model of the robot, either geometrical or of appearance [10,11], is learnt previously to the tracking task. In [12,13], the position of static and dynamic objects is obtained by multiple camera fusion inside an occupancy grid. An appearance model is used afterwards to ascertain which object is each robot. Despite the technique used for tracking, the common point of many of the proposals found in the topic comes from the fact that rich knowledge is obtained previously to the tracking, in a supervised task.
In this paper the localization of the robot is obtained in terms of its natural appearance. We propose a full metric and accurate system based on identifying natural features belonging to the geometric structure of the robot. On most natural objects we can find points whose image projection can be tracked in the image plane independently of the position the object occupies and based on local properties found in the image (i.e., lines, corners or color blobs). Those points are considered natural markers, as they serve as reference points in the image plane that can be easily relate with their three-dimensional counterparts. The set of methods focused on tracking natural markers have become a very successful and deeply studied topic in the literature [14,15], as they represent the basic measurements of most of existing reconstruction methods.
Scene reconstruction from image measurements is a classic and mature discipline in the computer vision field. Among the wide amount of proposals it can be highlighted those grouped under the name "Bundle Adjustment" [16,17]. Their aim is essentially to estimate, jointly and optimally, the 3D structure of a scene and the camera parameters from a set of images taken under some kind of motion. (i.e., it can be the camera that moves or equally some part of the scene w.r.t. the camera).
In general terms, Bundle Adjustment reconstruction methods are based in iterative optimization methods which try to minimize the image discrepancy between the measured positions of the 3D model and the expected ones using the last iteration solution. The discrepancy balances the contribution of the measurements into the final solution and plays an important role in this paper. Our main contribution is a redefinition of the discrepancy function using a Maximum Likelihood approach which takes into account the statistical distribution of the error. This distribution is especially affected by the odometry errors which are accumulative in long trajectories.
On the other hand, once a geometric model is obtained using a structure-from-motion approach, its pose with respect to a global coordinate origin can be easily retrieved by measuring the projection of the model in the image plane. This problem, commonly known as the Perspective n Point Problem (PnP), has received considerable attention in the literature, where some accurate solutions are found such as [18] or the recent global solution proposed in [19]. In this paper we instead follow a filtering approach, where not only image measurements but also last time pose information and odometry information are used to obtain the pose. This approach, which is based on the use of a Kalman Filter, is much more regular than solving the PnP problem at each time instant and will be described in this paper.
The paper is organized as follows. In Section 2. the notation and mathematical elements used in the paper are described. In Section 3. the description of the initialization algorithms is given. The online Kalman loop is explained in Section 4. and some results in a real environment are shown in Section 5. Conclusions are in Section 6.

Definitions and Notation
This section presents a description of the symbols and mathematical models used in the rest of the paper.

Robot and Image Measurements
The robot's pose at time k is described by a vector X k . We suppose that the robot's motion lie on the plane z = 0 referred from the world coordinate origin O W (See Figure 2). The pose vector X k is described by 3 components (x k , y k , α k ), corresponding to two position coordinates x k , y k and the rotation angle α k in the z axis. The motion model X k = g(X k−1 , U k ) defines the relationship between the current pose X k with respect to the previous time one X k−1 and the input U k given by odometry sensors inside the robot. In our proposal the robot used is a differential wheeled platform and thus U k = (V l k , Ω k ) T , where V l k and Ω k describe the linear and angular speed of the robot's center of rotation (O R in Figure 2) respectively. The motion model g is then very simple in function of the discretized linear speed V l k and the angular speed Ω k .
The robot's geometry is composed of a sparse set of N 3D points M = {M 1 , · · · , M N } referred from the local coordinate origin O R described by robot's pose X k . The points M i are static in time due to robot's rigidness, and thus, no temporal subindex is required for them. Function M i X k = t(X k , M i ) uses actual pose X k to express M i in the global coordinate origin O W that X k is referred to (see Figure 2): where: The augmented vector X a k , which is the state vector of the system, is defined as the concatenation in one column vector of both the pose X k and the set of static structure points M: An augmented motion model g a is defined as the transition of the state vector: It must be noticed that the motion model g a leaves the structure points contained in the state vector untouched, as we suppose that the object is rigid.
The whole set of measurements obtained in the image plane of the camera at time k is denoted by Y k .
Each measurement inside Y k is defined by a two-dimensional vector , describing the projection of each point inside M in the image plane, where i stands for the correspondence with the point M i . It must be remarked that, in general, Y k does not contain the image projection of every point in the set M as some of them can be occluded depending on the situation.
The camera is modeled with a zero-skew "pin-hole" model (see [17] for more details), with the following matrix K c containing the intrinsic parameters: The extrinsic parameters of the camera (i.e., position and orientation of the camera w.r.t. O W ) are described using the rigid transformation R c , T c (i.e., rotation matrix and translation vector). The matrix R c is defined by three Euler angles α c , β c , γ c . The vector T c can be decomposed into its three spatial coordinates T x , T y , T z . All calibration parameters can be grouped inside vector P : Given a single measurement y i k from the set Y k , its 2D coordinates can be expressed using the "pin-hole" model and the aforementioned calibration parameters: where λ i k represents a projective scale factor that can be obtained so that the third element of the right part of Equation (7) is equal to one. It is important to remark that the projection model, although simple, depends in a non-linear fashion w.r.t. M i due to the factor λ i k . For compactness the projection model of Equation (7) can be described as the function h: In the same way, the whole vector Y k can be expressed with the following function:

Random Processes
This paper explicitly deals with uncertainties by describing every process as a random variable with a Gaussian distribution whose covariance matrix stands for its uncertainty. The random processes are expressed in bold typography (i.e., X k is the random process describing the pose and X k is a realization of it) and each of them are defined by a mean vector and a covariance matrix. Therefore X a k is defined by its meanX a k and covariance Σ a k (or simplifying N (X a k , Σ a k )). The following processes are considered in the paper: , whose uncertainty comes from errors in image detection.
• Odometry input values U k = N (Û k , Σ U k ), which are polluted by random deviations.

Initialization Process
The initialization step consists of obtaining the initial pose the robot occupies and a 3D model of its structure using visual features (i.e., to obtain an initial guess of X a 0 ). The importance of this step is crucial, as the robot's geometric model serves as the necessary link between robot's pose and image measurements. The initialization allows to use afterwards the online approach presented in Section 4.
A delayed initialization approach is proposed, based on collecting image measurements and the corresponding odometry position estimation along a sequence of time (i.e., k = 1, · · · , K). After taking a sufficient amount of measurements of the object in motion, an iterative optimization method is used to obtain the best solution for X a 0 according to a cost criterion. The odometry readings are used in this step as metric information, which allows to remove the natural ambiguity produced by measurements from a single camera. The main problem of using odometry as a metric reference for reconstruction comes from the accumulative error it presents, which is a very well-known problem in dead-reckoning tasks. The initialization algorithm consists of a Maximum Likelihood (M.L.) estimation of the vector X a 0 , which does not estimate the initialization value itself but its equivalent Gaussian distribution X a 0 = N (X 1 0 , Σ a 0 ) by using a "Bundle Adjustment" method (See Figure 3). The M.L. approach (see [16,17] for more details) allows to properly tackle the growing drift present in the odometry estimation by giving fair less importance or weight to those instants where the odometry is expected to be more uncertain.
In the following sections, the initialization algorithm is explained in detail.

Odometry Estimation
A set of odometry readings, U 1 , · · · , U K are collected over a set of K consecutive time samples k = 1, · · · , K, which corresponds to the initialization trajectory performed by the robot. Using a motion model g given initial position X 0 , an estimation of the robot's pose in the whole initialization sequence is obtained as follows: . . .
where we recall that X k and U k denote Gaussian processes. Using expression (10), and by propagating the statistical processes through function g, the joint distribution p(X 1 , · · · , X K |X 0 ) can be obtained. The propagation of statistical processes through non-linear functions is in general a very complex task, as it requires to solve integral equations. However in this paper we will follow a first order approximation of the non-linear functions centered in the mean of the process (see [20] for more details). Using this technique the Gaussian processes X k and U k can be iteratively propagated through function g at the cost of being an approximation. This paper shows that despite being an approximation this approach end in a good estimation of the initialization vector. By denoting as X to the joint vector containing all poses X = (X T 1 , · · · , X T K ) T , the joint distribution of all poses, given the initial state X 0 , p(X|X 0 ) is approximated as Gaussian p.d.f. with meanX and covariance matrix Σ X .

Image Measurements
In this step we describe how to collect the position in the image plane of different points from the robot's structure during the initialization trajectory. We base our method in the work [15] where the SIFT descriptor is introduced. The process used in the initialization is composed of two steps: • Feature-based background subtraction.
We describe the static background of the scene using a single image, from which a set of features, The sets D b and F b are our feature-based background model.
Given an input image, namely I k , we find the sets D k and F k . We consider that a feature This method although simple shows to be very effective and robust in real environments (See Figure 4 ). Given the first image of the initialization trajectory, namely I 1 , we obtain the set of N features F 1 , once they are cleaned by the aforementioned background subtraction method. We propose a method to track them in the whole initialization trajectory.
By only using SIFT descriptors and its tracking in consecutive frames does not produce stable tracks, specially when dealing with highly irregular objects where many of the features are due to corners. We thus propose to use a classical feature tracking approach based on the Kanade-Lucas-Tomasi (KLT) algorithm [14]. To avoid degeneracy in the tracks, which is a very common problem in those kind of trackers, we use the SIFT descriptors to remove those segments of the tracks that clearly do not correspond to the feature tracked. This can be done easily by checking distance between the descriptors in the track. The threshold used must be chosen experimentally so that it does not eliminate useful parts of the tracks. In Figure 5a we can see the tracks obtained by the KLT tracker without degeneracy supervision. In Figure 5b the automatically removed segments are displayed.

Likelihood Function
Using the feature-based algorithm proposed before a set of measurements in the whole trajectory, Y 1 , · · · , Y K , is obtained, where each vector Y k contains the projection of N points from robot's structure at time sample k. The set of N points, M 1 , · · · , M N , jointly with the initial pose X 0 , represent the initialization unknowns. As the robot's motion does not depend of its structure, the following statistical independence is true: If all trajectories are packed into vector Y L , the following function put in relation Y L with the distribution of expression (11): where V k represents the uncertainty in image detection. Using distribution showed in (11), and propagating statistics through function (12), the following likelihood distribution is obtained: The likelihood function (13) is represented by a Gaussian distribution using a first order approximation of (12). It is thus defined by a meanŶ L and a covariance matrix Σ L :

Maximum Likelihood Approach
The likelihood function (14), parametrized by its covariance matrix Σ L and its meanŶ L , is dependent of the conditional parameters and unknown values X a 0 . The "Maximum Likelihood" estimation consists of finding the values for X a 0 that maximize the likelihood function: In its Gaussian approximation and by taking logarithms, it is converted into the following minimization problem: where Y L is the realization of the process (i.e., the set of measurements from the image) andŶ L are the expected measurements given a value of the parameters X 0 , M 1 , · · · , M N .
The configuration of Σ L is ruled by the expected uncertainty in the measurements and the statistical relation between them. Usually all cross-correlated terms of Σ L are non-zero, which has an important effect in the sparsity of the Hessian used inside the optimization algorithm and consequently in the computational complexity of the problem (see [21] for more details).
The covariance matrix Σ L can be approximated assuming either independence between time samples (discarding cross-correlation terms between time indexes) or total independence between time and each measurement (discarding all cross-correlation terms except the 2×2 boxes belonging to a single measurement). In Table 1, the different cost functions are shown depending on the discarded terms of Σ L . The different approximations of Σ L have a direct impact in the accuracy and the reconstruction error, and the results will be shown in Section 5.
Identity approximation of Σ L (M.I.) Intuitively, if Σ L results to be a identity matrix, the cost function is reduced to a simple image residual minimization extensively used in Bundle Adjustment techniques, where in principle all cost differences The result of minimizing Equation (16) instead of the usual (17) show significant improvements in the reconstruction error. In Section 5. a comparison is shown between the initialization accuracy under the different assumptions of the matrix Σ L , from its diagonal version to the complete correlated form. The minimum of (16) is obtained using the Levenberg-Mardquardt [17] iterative optimization method.
We suppose that during the measurement step there is a low probability of encountering outliers in the features. This argument can be very optimistic in some real configurations where more objects appear in the scene and produce occlusions or false matchings. For those cases, all cost functions presented in this section can be modified to be robust against outliers by using M-Estimators. We refer the reader to [17] for more details.

Initialization before Optimization
The use of an iterative optimization method for obtaining the minimum of (16) requires a setup value from which to start iterating. An initial estimation of X a 0 close to its correct value reduces the probability to fall in a local extrema of the cost function.
The method proposed in this paper to give an initial value for X a 0 consists of a non-iterative and exact method, which gives directly an accurate solution for X a 0 in absence of noise in odometry values. This method is based on the assumption that the robot moves in a plane and thus only the angle over its Z axis is needed. As explained below, this assumption allows to solve the problem with a rank deficient Linear Least Squares approximation, which is solved using the Singular Value Decomposition of the system's matrix and imposing a constraint that ensures a valid rotation. Its development is briefly introduced next.
For a point M i of robot's model and at time k of the initialization trajectory, the image measurement y i k results from the following projection: where the matrix R ∆ k and vector T ∆ k represent the rotation and position of the robot in the floor plane given by the odometry supposing that X 0 = ( 0 0 0 ) T . The rotation matrix R 0 and the offset T 0 correspond with the initial pose X 0 = (x 1 0 , x 2 0 , α 0 ) in form of rigid transformation: where R 0 is a non-linear function of the orientation angle. The expression (18) depends non-linearly of vector X a 0 and so a new parametrization is proposed jointly with a transformation which removes the projective parameter λ i k . The point M i is replaced by the rotated M i X 0 = R 0 M i , removing thus the product between unknowns. The orientation angle in the pose is substituted by parameters a = cos(α 0 ) and b = sin(α 0 ), with the constraint a 2 + b 2 = 1. Using the new parameterization the expression (18) is transformed in the following: with L T ∆ k : where x 1,∆ k y x 2,∆ k are the two first coordinates of T ∆ k . The new unknown vector, which correspond to the new parametrization of X a 0 , is: The expression (20) is decomposed in the following elements: where U i k , V i k y S i k are linear in terms of Φ but not in terms of y i k .
The projective scale is removed from the transformation by using vector product properties:  where two lineally independent equations are obtained for Φ.
Using all measurements inside Y L , a lineal system of equations is obtained in terms of Φ: It is straightforward to show that system of (27) has a single-parameter family of solutions. If Φ 0 is a possible solution, then Φ 0 + ψ∆ is a solution for any ψ ∈ R, with ∆: In fact, if ∆ is normalized, it matches up with the eigen-vector associated to the zero eigenvalue of matrix A T A.
Using the constraint that a 2 + b 2 = 1, and the singular value decomposition of matrix A, an exact solution of system (20) is obtained.
Once Φ is available, the solution X a 0 is obtained by inverting the parametrization: This method, although exact, is prone to error due to odometry noise and does not benefit from the Maximum Likelihood metric. However it is valid as a method to give an initial value for X a 0 before using the iterative approach.

Degenerate Configurations
The kind of trajectory performed by the robot during initialization has direct influence in the solution of X a 0 . There are three kinds of movements that yields degenerate solutions: • Straight motion: there is no information about the center of rotation of the robot and thus the pose has multiple solutions.
• Rotational motion around robot axis: the following one-parameter (i.e., n) family of solutions gives identical measurements in the image plane: • Circular trajectory: under a purely circular trajectory the following one-parameter family of initialization vectors gives identical measurements in the image plane: In practical cases it has been proved to be effective enough to combine straight trajectories with circular motion to avoid degeneracies in the solution of X a 0 .

Obtaining the Gaussian Equivalent of the Solution
Once the minimum of (16) is reached we suppose that the resulting value of X a 0 is the mean of the distribution X a 0 . This section describes how to also obtain the covariance matrix. The covariance matrix Σ a 0 of the optimized parameters is easily obtained by using a local approximation of the term Y L −Ŷ L in the vicinity of the minimumX a 0 using the following expression: where J is the Jacobian matrix ofŶ with respect to parameters X a 0 . The Jacobian is available from the optimization method, in which is used to compute the iteration steps.

Online Robot Localization
In this section the solution to X a k given the last pose information is derived. The fact that last frame information is available and the assumption of soft motion between frames allows to greatly simplify the problem.
A special emphasis is given to the fact that any process handled by the system is considered a random entity, in fact a Gaussian distribution defined at each case by its mean vector and covariance matrix. The problem of obtaining pose and structure, encoded in X a k given image observations Y k and the previous time estimation X a k−1 is viewed from the point of view of statistical inference, which means searching for the posterior probability distribution p(X a k |Y 1 , · · · , Y k ). That distribution gives the best estimation of X a k given all the past knowledge available. In Figure 6, a brief overview of the online method is presented. The online approach is divided into three steps: • Estimation Step: using the previous pose distribution p(X a k−1 |Y 1 , · · · , Y k−1 ), defined by its mean X a k−1 and covariance matrix Σ a k−1 , and the motion model g a , a Gaussian distribution which infers the next state is given p(X a k |Y 1 , · · · , Y k−1 ).
• Robust Layer: the correspondence between image measurements and the 3D model of the robot easily fails, so a number of wrongly correspondences or outliers pollute the measurement vector Y k . Using a robust algorithm and the information contained in the state vector, the outliers are discarded before the next step.

• Correction
Step: using an outlier-free measurement vector, we are confident to use all the information available to obtain the target posterior distribution p(X a k |Y 1 , · · · , Y k ) In all three steps it is required to propagate statistic processes over non-linear functions (f and h). As was stated in the previous section we show how to face the problem using first order expansions as it offers more compactness and is more readable. As a consequence the "Estimation" and "Correction" steps are solved using the so called Extended Kalman Filter (EKF) equations, which have been already implemented in problems of similar complexity such as visual Simultaneous Localization and Mapping (SLAM) [5].

Estimation Step
The estimation step uses the motion models available to infer the next pose of the robot which implies an increment in uncertainty.
Starting from the last pose distribution p(X a k−1 |Y 1 , · · · , Y k−1 ) = N (X a k−1 , Σ a k−1 ), the motion model g a and the noise included in odometry values, the following update is obtained: whereX a k|k−1 and Σ a k|k−1 are the mean and covariance matrix of distribution: The matrices J X and J U are the first derivatives of the function g a with respect to X a k−1 and U k respectively. Usually J X in odometry systems is the identity, therefore, at this step, the covariance matrix Σ a k|k−1 results to be bigger in terms of eigenvalues, which means uncertainty.

Correction Step
The correction step removes the added uncertainty in the estimation by using image measurements. It passes from the distribution p(X a k |Y 1 , · · · , Y k−1 ) to the target distribution p(X a k |Y 1 , · · · , Y k ), which includes the last measurement.
Using the estimation shown in (35), and knowing the correspondence between measurements with the camera and structure point of the state vector, the estimated measurement is given: where J h is the Jacobian matrix of the function h a with respect to X a k and Σ V is block diagonal matrix with Σ v on each block. Function h a performs the projection in the image plane of the camera of all visible points that form up the measurement vector Y k . The correction step itself is a linear correction of X a k|k−1 and Σ a k|k−1 by means of the Kalman gain K G : As it is stated in (41) the resulting Σ a k is reduced compared to Σ a k|k−1 which means that after the correction step, the uncertainty is "smaller".

Robust Layer
The robust layer has the objective of removing bad measurements from Y k to avoid inconsistent updates of X a k in the correction step. In this paper we propose to include the Random Sample Consensus (RANSAC) algorithm [22] between the estimation and correction step of the filter. The general idea is to found among the measured data Y k a set which agrees in the pose X k , using the 3D model obtained in X a k−1 . The interest of applying RANSAC in a sequential update approach resides on several reasons: firstly it allows to efficiently discard outliers from Y k preventing algorithm's degeneracy, which happens even if the motion model is accurate. Secondly, compared to online robust approaches, where a robust cost function is optimized, RANSAC allows not to break the Kalman filter approach, as it only cleans the measurement vector of outliers. Furthermore we have observed experimentally that the RANSAC algorithm can be very fast between iterations, as only a few outliers are inside the data. (We use the RANSAC implementation described in [17], which implements a dynamical computation of the outlier probability).
The RANSAC method proposed in the commented framework obtains the consensus pose X k from the set of measurements Y k and the 3D model available in X a k−1 using the algorithm presented in [18]. For a robot moving in a plane, as it is the case with the mobile robot considered in this paper, it is enough to use a minimum of 2 correspondences between the model and the measurement which makes the RANSAC very fast for removing outliers.

Image Measurements
Contrary to the initialization case, in this step we have an accurate prediction of the tracked points in the image plane at each time instant, namely vector Y k|k−1 . Using such prediction we can easily match the 3D points in the state vector with measurements taken in a measurement set F k , using the SIFT detector applied to current image I k . The matching is done in terms of distance in the image plane. Let y i k|k−1 a feature inside vector Y k|k−1 and y j k a candidate obtained using SIFT method. We conclude that they are matched if |y j k − y i k|k−1 | 2 M < τ max , where || M states for the Mahalanobis distance using the covariance Σ y k|k−1 computed from the matrix Σ Y k|k−1 . The Mahalanobis distance allows to take into account the uncertainty predicted for y i k|k−1 and also helps to select a threshold τ max with a probabilistic criterion.

Results
This section describes the experimental setup developed to support the theoretical algorithms proposed in this paper. The experiments are divided in those performed over synthetic data and those run in a real implementation of the Intelligent Space in the University of Alcala (ISPACE-UAH) [23]. In both kind of experiments the same camera parameters are used, derived from the real device used in its real placement. The single camera consists of a CCD based sensor with resolution of 640 × 480 pixels and a physical size of 1/2 (around 8 mm diagonal). The optical system is chosen with a focal length of 6.5 mm which gives around 45 o of Field of View (FOV). The camera is connected to a processing and acquisition node through a IEEE1394 port, which support 15 fps in RGB image format acquisition. The intrinsic parameters are the following: The camera is placed with respect a global coordinate origin, as it is displayed in Figure 7. Camera calibration is performed previously to the localization task, using checkerboards as calibration patterns and the "Matlab Calibration Toolbox" implemented by [24]. The distortion parameters of the camera are not considered in this case. As it can be shown in Figure 7, the camera is placed in oblique angle, which is specially useful for covering large areas with less FOV requirements (less distortion) compared to overhead configurations. The robotic platform is connected to the same processing node controlling the camera by means of a wireless network. The camera acquisition and the odometry values obtained from the robot are synchronized. The control loop of the robot is internal, and it is prepared to follow position landmarks based on odometry feedback at faster sampling frequency than the image localization system (15 fps). Therefore, for each image acquisition cycle, the odometry system obtain several readings that are used by the internal loop control. In this paper the localization information provided by the cameras is not included in the control loop of the robot.

Simulated Data
In this experiment the robot structure is simulated with a set of N = 10 points distributed randomly inside a cylinder with radius R = 0.5m and height h = 1m. The odometry system is supposed to have an uncertainty σ 2 V l = 10 and σ 2 Ω = 1 in linear and angular speed respectively. The initialization trajectory is shown in Figure 8. The image measurements are considered polluted with a Gaussian process of σ 2 v = 10, independently of each measurement and image coordinate. To compare the performance of the initialization method proposed in Section 3., the following error magnitudes are described: where T 0,gth , α 0,gth and M i gth correspond to the ground truth values of the initialization vector X a 0 . The following two experiments are proposed: so that it can be guessed the relationship between trajectory length and initialization errors. In Figure 10 the initialization errors are displayed versus ρ t .
In light of the results shown in Figure 10a-c, the M.C. method is capable of reducing the error no matter how large is the trajectory chosen. However, in the rest of the approximations of Σ L there is an optimal point where the initialization errors are minimum. This results make sense, as without statistical modelling large trajectories contain accumulative errors which usually affects the final solution of X a 0 . Both experiments clearly manifest that the complete matrix Σ L , with all its cross-correlated terms (M.C.), outperforms the rest of proposals, especially when it is compared to the case where Σ L is the identity matrix (M.I.), which means no statistical modeling.

Real Data
The initialization experiment proposed in this paper consists of a robot performing a short trajectory from which its 3D model and initial position is obtained using the results of Section 3.. We present a comparison of the initialization results using three different trajectories. Each one of the trajectories is displayed in Figure 11 as a colored interval of the whole trajectory performed by the robot in the experiment. The results of the initialization method on each of the 3 trajectories selected is shown in Figure 12. Depending on the trajectory used, the features viewed by the camera vary and thus the corresponding initialized 3D model. As it can be seen in Figure 12, the 3D model is accurate and its projection matches with points in the robot's structure in all cases. It must be remarked that on each case, the position obtained as a result of the initialization (i.e., X 0 ) corresponds to the first position of each interval.  The ground truth is obtained by means of a manually measured 3D model of the robot. The model is composed of 3D points, that are easily recognized by a human observer. By manually clicking points of the 3D model in the image plane, and by using the method proposed in [18], the reference pose of the robot is obtained in some of the frames of the experiment.
Another experiment is proposed to test the online algorithm with occlusions and a larger path followed by the robot. In Figure 14a it is shown the geometric placement of the camera and in Figure 14b it is shown the geometric model obtained during initialization.  The resulting localization error is shown in Figure 15b. In Figure 16 it is shown several frames, indexed by the time sample number k, presenting hard occlusions between the camera and the robot without losing the tracking. The RANSAC method used in the Kalman loop avoid erroneous matches in the occluded parts of the object. Besides, in those frames with completely occluded features, only the estimation step of the Kalman filter is done, which is accurate enough for short periods of time. In both sequential experiments the 3D model is composed of 8 to 10 points and no extra geometrical information is introduced in the state vector. As a future line of study a simultaneous reconstruction and localization approach can be adopted to introduce online extra 3D points in the state vector as the robot is tracked. A very similar approach is done in the Simultaneous Localization and Mapping (SLAM) problem, and some of their solutions [5] are perfectly applicable to the problem assessed in this paper.

Conclusions
This paper has presented a complete localization system for mobile robots inside intelligent environments with a single external camera. The objectives of obtaining the pose of the robot based on its natural appearance has been tackled in the paper using a reconstruction approach followed by a sequential localization approach.
The main contributions of this paper are summarized in the following list: • The initialization step provides a non-supervised method for obtaining the initial pose and structure of the robot previously to its sequential localization. A Maximum Likelihood cost function is proposed, which obtains the pose and geometry given a trajectory performed by the robot. The proposal of this paper allows to compensate for the odometry drift in the solution. Also an exact initialization method has been proposed and the degenerate configurations have been identified theoretically.
• The online approach of the algorithm obtains the robot's pose by using a sequential Bayesian inference approach. A robust step, based on the RANSAC algorithm, is proposed to clean the measurement vector out of outliers.
• The results show that the proposed method is suitable to be used in real environments. The accuracy and non-drifting nature of the pose estimation have been also evaluated in a real environment.
The future research must be oriented to scale the problem for large configurations of non-overlapped cameras and multiple robots, where extra problems arise, such as to automatically detect the transitions between cameras and online refinement of the geometric models. It is important from our point of view to, in the future, combine the information given by the system proposed in this paper with information sensed onboard the robots using cameras. This approach allows to jointly build large maps attached to information given by the cameras, so that robots can be localized and controlled to perform a complex task.