convenience systems, safety systems, productivity systems and trac-assist systems, which are listed as follows:

• Convenience Systems Parking Assistance Adaptive Cruise Control (ACC) Lane-Keeping Assistance Automated Vehicle Control http://telematicsnews.info/2013/01/16/naias-mercedes-benz

During the past decades, due to the reason of road safety, advanced driver assistance systems (ADAS) or intelligent vehicles (IV) have obtained more and more attentions and developments from research society and industry community. Advanced driver assistance systems (ADAS), are systems to help driver in the driving tasks. ADAS usually consists of adaptive cruise control (ACC), lane departure warning system, collision avoidance system, automatic parking, trac sign recognition, etc. ADAS provides relatively basic control assistance by sensing the environment or assessing risks. Nowadays, ADAS have already been applied in some top-brand vehicles, e.g. Mercedes-Benz E-class . Intelligent vehicle (IV) systems are seen as a next generation beyond current ADAS. IV systems sense the driving environment and provide information or vehicle control to assist the driver in optimal vehicle operation. The range of applications for IV systems are quite broad and are applied to all types of road vehicles cars, heavy trucks, and transit buses. In [Bishop 2005], IV applications can generally

• Safety Systems

Assisting driver perception Crash prevention Degraded driving Precrash In intelligent vehicles, perception systems are of the most importance. Perception systems could sense and interpret surrounding environment based on various kinds of sensors, such as radar, 2D/3D lidar (laser range nder), monocular/binocular/omnidirectional vision systems, etc. Indeed, as described in [Fletcher 2008[START_REF] Petrovskaya | [END_REF]], all the experimental vehicles joined in 2007 DARPA grand challenge were equipped with advanced sensors (shown in Fig. 1.1).

In most cases, these sensors shoulder heavy responsibilities of environmental perception (self-localization, collision avoidance, motion planning, etc.). Towards the objective of environmental perception, stereo vision system and lidar are two conventional sensors. In general, stereo vision system is a passive sensor system consisting of at least two cameras. It oers abundant texture and potential semantic information of surrounding environments, and it is superior in 3D measuring at a relative low expense. Many perceptional tasks could be solved by stereo vision system, e.g. visual SLAM [START_REF] Davision | [END_REF]], visual detection and tracking [Srinivasa 2002], visual classication and recognition [Gavrila 2004]. On the other side, lidar is an active range sensor which eciently measures surrounding distances by visible/invisible laser. The eective detection range of lidar is always larger than stereo vision system. Merely based on lidars, a great number of techniques are also developed for scene perception and understanding, e.g. [Streller 2002[START_REF][END_REF]. In this thesis, we focus on the usages of stereo vision system and lidar.

This thesis is part of the project CPER "Intelligence du Véhicule Terrestre" (Intelligence of ground vehicle), developed within Systems and Transportation Laboratory (SeT) of Institute of Research on Transportation, Energy and Society (IRTES), UTBM, France. The problem addressed in this thesis is how to map dynamic environment by stereo vision and lidar. The main objectives and contributions are as follows:

• The rst objective is to estimate the ego-motion of a vehicle when it is driven.

Ego-motion estimation is a fundamental part for an intelligent vehicle. Only by knowing the ego-motion of the vehicle itself at rst, an intelligent vehicle can further analyze the surrounding dynamic environment. To achieve this objective, we implement a stereo vision based ego-motion estimation method.

Our contribution is a comparison between tracker based and matching based circular feature associations as well as performance analysis of image feature detectors. From the comparison, the best image feature detector and feature association approach are selected for estimating ego-motion of the vehicle.

• The second objective is to detect and recognize independent moving objects around the moving vehicle. This work is important for collision avoidance. We propose a stereo vision based independent moving object detection method, as well as a spatial information based object recognition method. The independent moving objects are detected based on the results of ego-motion estimation, with a help of U-disparity map. The recognition is performed by training classiers from spatial information of the moving objects.

• The third objective is extrinsic calibration between a stereoscopic system and 2D lidar. The extrinsic calibration between the two sensors is to calculate a rigid 3D transformation, which connects the positions of the two sensors. We use an ordinary chessboard to achieve the extrinsic calibration and intrinsic calibration of the stereoscopic system at the same time. Also, we present an optimal extrinsic calibration of the stereoscopic system and the lidar basing on geometric constraints that use 3D planes reconstructed from stereo vision.

• The fourth objective is mapping the sensed environment by occupancy grids. Occupancy grid map is a classic but practical tool of interpreting environment. In the thesis, we propose a stereo vision based occupancy grid mapping method. Our contribution is about the improvement of mapping results by analyzing the ground plane and augmenting the mapping by integrating motion and object recognition information. Then, we attempt to combine stereo vision and lidar to improve the occupancy grid map. As demonstrated in Fig. 1.2, our platform SetCar is equipped by many sensors (2D lidar, stereo vision system, sheye camera, RTK-GPS, gyro, etc). A Bumblebee XB3 stereo vision system is installed on the top, and oriented to the front. It contains three collinear cameras with a maximum baseline of 240mm. Images are captured in a format of 1280 × 960 pixels and a speed of 15 fps (frames per second). In our application, only the left and right cameras are used. A SICK LMS 221/291 lidar provides range measurements in a 180 • scanning plane, with an angular resolution of 0.5 • . It scans in 75Hz and 80 meters as maximum detective range. All of the installed sensors are xed rigidly. 1.4 Structure of the manuscript

In chapter 2, we review basic knowledges used in the thesis. The basic knowledges contain the measurement models of stereo vision system and lidar, linear/nonlinear least square methods and image feature detectors/descriptors. The reviewed basic knowledges make understanding the thesis more easy.

In chapter 3, we rst summarize a framework of stereo visual odometry. Within the framework, dierent kinds of feature detection and association are then compared at length. From the comparison results, the best approaches of detecting and associating image features are selected.

In chapter 4, independent moving objects are rstly detected and segmented based on the results of ego-motion estimation described in chapter 3. Next, spatial information of the segmented objects are extracted based on a kernel PCA method. Several classiers are learned from these spatial information for achieving independent moving object recognition.

In chapter 5, we propose a chessboard based extrinsic calibration between a stereoscopic system and a 2D lidar. To improve the accuracy of calibration results, we model the sensor measurements and take the sensor models into account. We show that by considering the sensor models, the calibration results are improved compared with a previous method.

In chapter 6, dynamic occupancy grid map is built rstly by stereo vision system. It integrates motion and object recognition information. The mapping process relies on 3D reconstruction of the stereo measures and is improved by analyzing geometrical feature of ground plane. At the same time, occupancy grid map is also created by lidar measures. By the extrinsic calibration results in chapter 5, we adopt a linear opinion pool to combine the two kinds of measures to create a occupancy grid map.

In chapter 7, conclusions and some research perspectives for this thesis are presented.

+D=FJAH

Basic Knowledge As demonstrated in Fig. 1.2, our platform SetCar is equipped by many sensors (2D lidar, stereo vision system, sheye camera, RTK-GPS, gyro, etc). In this section, we introduce the sensor measuring models of lidar and stereo vision system repectively.

Coordinate Systems

To analyze such a multi-sensor system, we set several coordinate systems with respect to the dierent sensors, as demonstrated in Fig. 2.1. Let R 3 stereo , R 3 lef t , R 3 right , and R 3 lidar , be the coordinate systems attached to the stereo vision system, the left and right cameras, and the lidar respectively. All the coordinate systems are connected with each other by a rigid Euclidean transformation (a rotation matrix R and a translation vector T ). For the reason of simplicity, we assume the left camera coordinate system as the reference of the stereo vision coordinate system, hence, R 3 stereo = R 3 lef t . After setting up all the sensor coordinate systems and their connections, we are going to introduce measurement models of these sensors. In the lidar coordinate system, we specify the origin as the point which emits laser rays. The directions X, Y , Z are set as rightward, upward and forward from the sensor respectively, as shown in Fig. 2.1. The lidar we used (SICK LMS 221/291) provides an array of raw distance measurements r i , θ i in polar coordinate system. In real applications, these raw measurements are converted into Cartesian form (X i Y i Z i ) T (note that in our assumption, Y i always equals to zero):

Stereo Vision System

[ X i Z i ] = [ r i × cos θ i r i × sin θ i ] (2.1)
Error Model: In most cases, lidar measurements contain errors, as described in the manual of SICK LMS 221/291 . According to this manual, we assume the lidar measurements, ranges r i and scan angles θ i , have additive Gaussian noises:

r i = R i + n r θ i = Θ i + n θ (2.2)
where R i and Θ i are the "true" measurements. n r and n θ are independent additive zero-mean Gaussian noises with variances δ r and δ θ , respectively. Hence, after being converted into Cartesian form, we have: 2.3) http://sicktoolbox.sourceforge.net/docs/sick-lms-technical-description.pdf where Pi lidar is the "true" measurement in 3D Cartesian coordinate, n p is the variance vector. When assuming n θ ≪ 1 (which is the fact in most lidar equipments), hence, sin n θ ∼ n θ , cosn θ ∼ 1. Eq. ( 2.3) yields:

[ X i Z i ] = (R i + n r ) [ cos(Θ i + n θ ) sin(Θ i + n θ ) ] = (R i + n r ) [ cos Θ i cos n θ -sin Θ i sin n θ sin Θ i cos n θ + cos Θ i sin n θ ] = Pi lidar + n p (
n p = (R i + n r )n θ [ -sin Θ i cos Θ i ] + n r [ cos Θ i sin Θ i ] (2.4)
From the above, the covariance matrix for a LIDAR measurement is:

Q i lidar ≜ E[n p n T p ] = (R i ) 2 δ 2 θ 2 [ 2 sin 2 Θ i -sin 2Θ i -sin 2Θ i 2 cos 2 Θ i ] + δ 2 r 2 [ 2 cos 2 Θ i sin 2Θ i sin 2Θ i 2sin 2 Θ i ] (2.5)
where Q i lidar is the covariance matrix of the ith lidar measurement P i lidar . δ r and δ θ are the variances of independent additive zero-mean Gaussian noises n r and n θ . In practice, Θ i and R i are approximated by θ i and r i , respectively. Thus, in homogeneous coordinates, the covariance matrix of Pi lidar can be expressed as:

Σ Pi lidar =     Q i 11 0 Q i 12 0 0 0 Q i 21 0 Q i 22 0     (2.6)
where Q i kj corresponds to the (k, j)th element in Eq. (2.5). The error transformation is shown in Fig. 2 The imaging process for an ordinary CCD/CMOS camera is a transformation from photons into electrons. For most cameras with CCD/CMOS like sensors, pinhole model is a basic but useful model to present this mapping. A typical pinhole camera model is shown in Fig. 2.3. Let camera center O be the origin of camera coordinate system R 3 camera , and consider the image coordinate system R 2 image , which locates in the image plane Z = f , where f is focal length, the distance between the camera (2.7) Eq. 2.7 is inconvenient in practice. It can be written in terms of matrix manipulation under homogeneous coordinates:

Z   x y 1   =   f 0 0 0 0 f 0 0 0 0 1 0       X Y Z 1     (2.8)
Note that in homogeneous coordinates , a point (x, y, w) T is equivalent to (x/w, y/w, 1) T , where w is a scale factor. Hence,

  x y 1   =   f 0 0 0 0 f 0 0 0 0 1 0       X Y Z 1    
(2.9)

In fact, the coordinates in the image do not correspond to the physical coordinates in the image plane. For a CCD camera, the relationship depends on the size and shape of the pixels and of the position of the CCD chip in the camera. The transformation is illustrated in Fig. 2.4: The image coordinates (u, v) are obtained by: https://en.wikipedia.org/wiki/Homogeneous_coordinates Figure 2.4: From physical coordinates to image coordinates in image plane

  u v 1   =   f k u (tan θ)f k v c u 0 0 f k v c v 0 0 0 1 0       X Y Z 1    
(2.10)

where k u and k v are scale factors with [pixels/length] units, which are the number of pixels within a unit length along u and v direction, respectively. (c u , c v ) is the position of principal point. θ is the skew angle between u and v axis, as indicated in Fig. 2.4. Let f u = f k u , f v = f k v , and the skew angle θ = 0, Eq. 2.10 is simplied as:

  u v 1   =   f u 0 c u 0 0 f v c v 0 0 0 1 0       X Y Z 1    
(2.11)

Now, let

K =   f u 0 c u 0 f v c v 0 0 1  
(2.12) then, Eq. 2.11 is denoted as: 2.13) where p and P are the homogeneous coordinates of a point in image coordinate system and camera coordinate system, respectively. K is called intrinsic matrix, or camera calibration matrix. I is a 3 × 3 identity matrix.

p = K[I|0] P (

Lens Distortion Model

The pinhole camera model described previously in Sec. 2.1.3.1 is an ideal model without considering lens distortions. As a result of several types of imperfections in design and assembly of lenses composing the camera optical system, the expression of Eq. 2.13 does not hold true. In this section, we discuss about two types of lens distortion: radial distortion and tangential distortion. • Tangential distortion: This type of distortion is caused by errors in optical component shapes and optics assembly. It is also called decentering distortion.

It is illustrated in Fig. 2.5 (b) as dt.

Let a point (x d , y d ) be the distorted coordinates of a point (x, y) in image plane coordinate system R 2 image , and r 2 = x 2 + y 2 . Then:

( x d y d ) = (1 + k 1 r 2 + k 2 r 4 + k 3 r 6 ) • ( x y
) + dt (2.14) where k 1 , k 2 , k 3 are the coecients of radial distortion. dt is the tangential distortion vector: .15) where p 1 , p 2 are the coecients of tangential distortion. Therefore, under the lens distortion model, the pixel coordinates (u, v) are related to distorted coordinates similar to Eq. 2.13: 2.16) where (x d , y d ) is represented in Eq. 2.14. model, with intrinsic matrices K l and K r respectively. Therefore, suppose there is a point P in 3D space, its corresponding coordinates in the left camera coordinate system R 3 lef t , the right camera coordinate system R 3 right , the left image coordinate system R 2 lef t and the right image coordinate system R 2 right are P l , P r , p l , p r respectively. Then, they are connected as:

dt = [ 2p 1 xy + p 2 (r 2 + 2x 2 ) p 1 (r 2 + 2y 2 ) + 2p 2 xy ] ( 2
  u v 1   = K   x d y d 1   (
pl = K l [I|0] • Pl pr = K r [I|0] • Pr P l = R l r • P r + T l r (2.17)
where • denotes homogeneous coordinates. R l r , (a 3 × 3 rotation matrix) together with T r r (a 3×1 translation vector) represent rigid 3D Euclidean transformation from

R 3
right to R 3 lef t . In similar, R r l and T r l describe a rigid 3D Euclidean transformation from R 3 lef t to R 3 right . Usually, R and T are called extrinsic parameters. Another characteristic of binocular vision system is the epipolar geometry. In epipolar geometry, the point intersections of the base line O l O r with the two image planes are epipoles. Any point P in 3D space with two camera centers O l and O r dene an epipolar plane. The line intersections of the epipolar plane with the two image planes are called epipolar lines. Indeed, all epipolar lines intersect at epipole. An epipolar plane denes the correspondence between the epipolar lines. All the epipoles, epipolar plane and epipolar lines are drawn in Fig. 2.6. Based on the above denitions, epipolar constraint is that: suppose p l is the left image position for a scene point P , then, the corresponding point p r in the right image must lie on the epipolar line. The mathematical interpretation is that, for a corresponding point pair (p l , p r ) in homogeneous coordinates •, we have:

pT l F pr = 0 (2.18)
where F is called fundamental matrix. F pr describes a line (an epipolar line) on which the corresponding point pl must lie. Intrinsic calibration of a camera has been researched for decades. Among all the proposed methods, 2D planer pattern based methods, e.g. Zhengyou Zhang's method [Zhang 2000], is the most popular. This method requires the camera to observe a planar pattern (e.g. 2D chessboard) shown at a few dierent orientations. Either the camera or the planar pattern can be freely moved. It is implemented in a well known camera calibration toolbox for Matlab ! , developed by Jean-Yves Bouguet [Bouguet 1999]. This Matlab calibration toolbox is applied in our practices. The whole process is as follows:

• Prepare a chessboard pattern and attach it to a planar surface.

• Take a few images of the chessboard under dierent positions with dierent orientations, as shown in Fig. 2.7 (a).

• Detect the corner points on the chessboard images, as shown in Fig. 2.7 (b).

• Estimate the intrinsic parameters and the coecients of radial distortion through a nonlinear optimization process. The calibration results and error analysis are shown in Fig. 2.7.

Extrinsic Calibration of Binocular Vision System

Extrinsic calibration of a stereo rig is to estimate extrinsic parameters between left and right cameras (a rotation matrix and a translation vector, as described in Sec. 2.1.4). Given corresponding image points from two views, the extrinsic parameters can be calculated from SVD factorizing the fundamental matrix F [Hartley 2004]. Whereas, in Jean-Yves Bouguet's Matlab calibration toolbox, intrinsic parameters for the left and right cameras are calculated at rst. Then, the method in [Hartley 2004] is adopted for acquiring an initial estimation of the extrinsic parameters. At last, this initial estimation is optimized by minimizing reprojection errors of all the corresponding corner points in the two views. A demonstrative example in Fig. 2.8 shows the extrinsic calibration results of our stereo vision system. Image undistortion is to compensate the eects of lens distortions described in Sec. 2.1.3.2. Indeed, it consists of a transformation from physically distorted image to an ideal image under the pinhole camera model (see in Sec. 2.1.3.1). For every pixel in undistorted image, we have to compute its distorted location. That is, for each pixel (u, v) in the corrected image, undistortion process computes its corresponding x

← (u -c x )/f x y ← (v -c y )/f y x ′ ← x(1 + k 1 r 2 + k 2 r 4 + k 3 r 6 ) + 2p 1 xy + p 2 (r 2 + 2x 2 ) y ′ ← y(1 + k 1 r 2 + k 2 r 4 + k 3 r 6 ) + p 1 (r 2 + 2y 2 ) + 2p 2 xy map x (u, v) ← x ′ f x + c x map y (u, v) ← y ′ f y + c y (2.19)
where (c x , c y ), f x , f y are the intrinsic parameters calculated in Sec. 2.2.1.

k 1 , k 2 , k 3 , p

Stereo Rectication

Stereo rectication is to align the image planes of the two views. In practice, if the two cameras are aligned to be coplanar, matching point between two images is simplied to searching in one dimension -a horizontal line parallel to the baseline between the cameras. In fact, stereo rectication is to move the epipoles to innity and match up epipolar lines. Its eects are demonstrated in Fig. 2.9. An outline of stereo rectication is given in [Hartley 2004]:

1. Find initial point correspondences For the inputs of a stereo vision system (the image pairs acquired by the left and right cameras), image undistortion and stereo rectication are usually utilized as preprocessings to simplify 3D computer vision tasks. Be similar to Eq. 2.19, for each pixel (u, v) in the destination image, undistortion together with stereo rectication compute its corresponding coordinates in the original image:

x ← (u -c ′ x )/f ′ x y ← (v -c ′ y )/f ′ y [XY W ] T ← R -1 × [xy1] T x ′ ← X/W y ′ ← Y /W x ′′ ← x ′ (1 + k 1 r 2 + k 2 r 4 + k 3 r 6 ) + 2p 1 x ′ y ′ + p 2 (r 2 + 2x ′2 ) y ′′ ← y ′ (1 + k 1 r 2 + k 2 r 4 + k 3 r 6 ) + p 1 (r 2 + 2y ′2 ) + 2p 2 x ′ y ′ map x (u, v) ← x ′′ f x + c x map y (u, v) ← y ′′ f y + c y (2.20)
where (k 1 , k 2 , k 3 , p 1 , p 2 ) are the distortion coecients in one view (e.g. the left camera). R is the rectication transformation (3 × 3 matrix) for one view. c x , c y , f x , f y are the original intrinsic parameters for one view. c ′

x , c ′ y , f ′ x , f ′ y are intrinsic parameters of new camera projection matrix after rectication. Note that after rectication, the intrinsic matrices K l and K r are the same. That is, after undistortion and rectication, each stereo image pair can be thought as acquired by one camera with pure translation movement under ideal pinhole model.

The results after undistortion and rectication are shown in Fig. 2.10, where the left and right images are well aligned (the greelines are helped to see whether the corresponding points in the two views are colinear).

2.2.4

Corner Points Triangulation 3D triangulation is a process of reconstructing a point in 3D coordinates from its corresponding projections in two or more images, giving the conguration parameters of the cameras. Here, we only discuss 3D point triangulation from two views. .11: 3D triangulation between two parallel, identical cameras are assumed to be identical and the images planes are colinear. Let f be the focal length, b be the distance between the two camera centers, XZ is a plane where the optical axes lie, XZ plane is perpendicular to the image plane of both cameras, X axis equals the baseline. A 3D point P 's projection in XZ plane is P xoz = (X, Z). As marked in the gure, the two lines that connect connected camera centers and P intersect image planes at x l and x r , which are indeed the image coordinates in X axis of P 's two projections. It is easy to deduct the followings from similar triangles:

Z f = X x l Z f = X -b x r (2.21)
As the Y -axis is perpendicular to the XZ plane, we can have the following equation in similar:

Z f = Y y l = Y y r (2.22)
Usually, we name ∆ = x lx r as disparity. Then, 3D coordinates (X, Y, Z) of a point P in the left camera coordinate system can be derived as:

X = (u l -c u ) • b ∆ Y = (v l -c v ) • b ∆ Z = f • b ∆ (2.23)
where u l and v l are the pixel coordinates of P in the left image, (c u , c v ) is the position of left camera's principal point.

Least Squares Estimation

The method of least squares [START_REF] Rao | [END_REF]] is about estimating parameters of a model by minimizing the squared discrepancies between observed data and their expected values provided by a model. Suppose a data set consisting of n data pairs (x i , y i ), i = 1, ..., n, where x i is an independent observed variable vector and y i is an observed response whose value is decided by observation. The model is represented as

y = f (x, β) (2.24)
where β is a vector of parameters dened by the model. The goal is to nd β for the model which "best" t the observation. The least squares method aims to nd the optimal values of β by minimizing squared residuals:

S = n ∑ i=1 r 2 i r i = y i -f (x i , β) (2.25)
where S is a summation of residuals r i . According to whether or not the residuals r i are linear or not, least squares estimation can be divided into two categories: linear least squares and non-linear least squares. If measurement errors exist in observed independent variables x i , errors-in-variables [START_REF] Rao | [END_REF]] model may be considered instead of standard least square estimation. In this section, linear least squares/non-linear squares estimators and their variations are discussed.

Linear Least Squares

In linear least square problems, the model f (x i , β) can be expressed linearly in terms of parameters β. Mathematically, linear least squares solves an over-determined system of linear equation. Usually, linear least squares problem have a closed-form solution. For example, consider an over-determined system:

n ∑ j=1 X ij β j = y i , (i = 1, 2, ..., m) (2.26) of m linear equations with n coecients, β 1 , β 2 , ..., β n (m > n).
This can be written in matrix form as:

Xβ = y
(2.27) where

X =      X 11 X 12 . . . X 1n X 21 X 22 . . . X 2n . . . . . . . . . . . . X m1 X m2 . . . X mn      , β =      β 1 β 2 . . . β n      , y =      y 1 y 2 . . . y m     
(2.28)

As mentioned before, least squares estimation seeks the coecients β minimizing the summation of residuals:

β = arg min β S(β) S(β) = ∥y -Xβ∥ 2
(2.29)

The solution for Eq. 2.29 can be proven to be [START_REF] Rao | [END_REF]]:

β = (X T X) -1 X T y = X + y
(2.30) where X + is the Moore-Penrose pseudo-inverse [Bretscher 1995] of X.

Ordinary Linear Least Squares

When applying least squares in real data sets, errors/noises are inevitable. Suppose there exist errors between the actually observed responses y i and the "predicted" outcomes x T i β, the model in Eq. 2.27 can be written as:

y = Xβ + ε
(2.31) where y and ε are m × 1 vectors representing actual response values and errors. Assume the errors are zero mean independent and identically distributed (i.i.d) Gaussian variables that: ε ∼ N (0, σ 2 I), where I is an m × m identity matrix, and σ 2 determines the variance of each observation. The ordinary least squares estimation of β can be given by explicit formula as in Eq. 2.30:

β = (X T X) -1 X T y (2.32)
Hence, the predicted value from the observation is:

ŷ = X β = X(X T X) -1 X T y (2.33)
Beside estimating the value of β, the covariance matrix of least squares estimator is also important. It can be shown that, given X, the covariance matrix Q β of the estimator β, equals to [Hayashi 2011]:

Q β = (X T X) -1 σ 2
(2.34) σ 2 can be estimated from the data as:

σ2 = 1 m -n ∥y -X β∥ 2 (2.35)
Therefore, Q β can be estimated by:

Q β = (X T X) -1 σ2
(2.36) The ordinary linear least squares assumes that each data point provides equally precise information about the estimation. In other words, the error term ε is constant over all observations. Clearly, this assumption does not hold in most modeling applications. Weighted linear least squares oers a better approach by attempting to give each data point its proper amount of inuence over the parameter estimation.

In this case, one can minimize the weighted sum of squares:

β = arg min β m ∑ i=1 w i |y i -x i β| 2 = arg min β ||W 1/2 (y -Xβ)|| 2 = arg min β (y -Xβ) T W(y -Xβ) (2.37)
where w i > 0 is the weight of the ith observation, and W is a m×m diagonal weight matrix. Ideally, the weight matrix should be equal to the reciprocal of the variance of the measurement.

W = cov(y) -1 =      w 1 w 2 . . . w m      =       1 σ 2 1 1 σ 2 2 . . . 1 σ 2 m       (2.38)
where σ 2 i is the variance of error in ith observation data point. Indeed, W is the covariance matrix of the error ε. When W = 1/σ 2 • I, the weighted linear least squares becomes ordinary linear least squares. A closed-form solution is [Strutz 2010]:

β = (X T WX) -1 X T Wy (2.39)
The covariance matrix Q β of estimation β, can be acquired by applying error propagation law:

Q β = (X T WX) -1 X T WX(X T WX) -1
(2.40) However, in many real practices, we do not know σ 2 i . We could set the weights manually or by an iterative approach: iteratively reweighted least squares (IRLS) [START_REF] Chartrand | [END_REF]]. IRLS uses an iterative approach to solve a weighted least squares problem of the form:

β(t+1) = arg min β m ∑ i=1 w i |y i -x i β i | 2 (2.41)
In each step t + 1, the IRLS estimates β(t+1) by the aforementioned weighted linear least squares:

β(t+1) = (X T W (t) X) -1 X T W (t) y
(2.42) where W (t) is a diagonal weight matrix in tth step with each element:

w (t) i = |y i -x i β (t) | -1 (2.43) 2.3.1.3
Total Linear Least Squares

The above ordinary and weighted linear least squares only assume errors in the observed response y, as interpreted in Eq. 2.31. However, it is a more natural way to model data if both observation X and response y are contaminated by errors. Total least squares (TLS) [Huel 1991], is a natural generalization of ordinary least squares where errors in all observational variables are taken into account. TLS problem is dened as follows: we are given an over-determined set of m equations x i β = y i (i = 1...m) in n unknowns β j (j = 1...n), compiled to a matrix equation Xβ ≈ y. Both the vector y as well as the matrix X are subject to errors. The total least squares problem consists in nding an m × n matrix X and an m × 1 vector ŷ for which the equation

Xβ = ŷ
(2.44) has an exact solution, under the condition that the deviation between (X|y) and ( X|ŷ) is minimal in terms of the Frobenius norm:

min ∥(X|y) -( X|ŷ)∥ F (2.45)
For a m × n matrix A = (a i,j ) i = 1...m, j = 1...n, the Frobenius norm is dened as:

∥A∥ F = m ∑ i=1 n ∑ j=1 |a ij | 2 = √ trace(A * A) = min(m,n) ∑ i=1 σ 2 i (2.46)
where A * denotes the conjugate transpose, σ i are the singular values of A. Furthermore, we denote the errors lying in X and y as E (m × n matrix) and r (m × 1 vector), respectively. The TLS model becomes:

(X + E)β = (y + r) (2.
E,r ∥[E|r]∥ F subject to (X + E)β = y + r (2.48)
where [E|r] is an augmented matrix and ∥ • ∥ F is the Frobenius norm.

Let C = [X|y] be the m×(n+1) augmented matrix of X and y side by side, and C = UΣV T is its SVD decomposition. Assume that the singular values of C are such that: σ 1 ≥ ... ≥ σ k > σ k+1 = ...σ n+1 . Then, it has been shown in [Huel 1991] that β T LS minimizing Eq. 2.48 when all the errors are i.i.d Gaussian variables, is:

β T LS = - y α (2.49)
where the vector (y, α) of norm 1, with α ̸ = 0, is in the subspace S k spanned by the right singular vectors v k+1 , ..., v n+1 of V. To put it simple, β T LS is given by taking the negative of the results of normalizing the rst n elements of the last column V by the last element of V.

The ordinary least squares and total least squares methods evaluate the estimation accuracy in dierent ways: the ordinary least squares minimizes the sum of distances from data to the tting line, while the total least squares minimizes the sum of orthogonal distances instead. The dierence is shown in Fig. 2.12 However, in [Huel 1991[START_REF] Markovsky | [END_REF], the covariance of TLS estimator is not given. One important reason is that the likelihood function of TLS problem contains the true (noiseless) values of measurements as nuisance parameters [START_REF] Gleser | [END_REF]].

As a result, the dimension of likelihood function increases rapidly with the number of observations. In [START_REF] Nestares | [END_REF]], the likelihood function of TLS is simplied with several assumptions and the Cramer-Rao lower bound (CRLB) " of TLS is given:

Q (β) T LS ≥ CRLB(β T LS ) = 1 γ σ 2 n ∥(β) T LS ∥ 2 (X T X) -1 (2.50)
where σ n is the variance of the assuming i.i.d distributed measurement noises. γ =

σ 2 0 σ 2 n +σ 2 0
is related with the signal to noise ratio (SNR). If the SNR is high (σ

2 0 ≫ σ 2 n ), then γ ≈ 1.
" http://en.wikipedia.org/wiki/Cramer-Rao_bound

Non-linear Least Squares

Be dierent to the linear least squares introduced in Sec. 2.3.1, the observation model f (x, β) (in Eq. 2.24) in non-linear least squares can not be represented in form of linear combination of parameters β = β i (i = 1...n). Solutions of nonlinear least squares usually approximate non-linear model by a linear one, and rene the parameter vector β iteratively. Recall to Eq. 2.24, given m data points (x 1 , y 1 ), (x 2 , y 2 ), ..., (x m , y m ) and a non-linear model y = f (x, β), with m ⩾ n. Non-linear least squares is to nd the parameter vector β, which minimizes the sum of squares:

S = m ∑ i=1 r 2 i (2.51)
where the residuals r i are given by:

r i = y i -f (x i , β), i = 1, 2, ..., m (2.52)
According to the knowledge of calculus, the sum of squares S reaches a minimum when the gradient is zero. Since the model contains n parameters, there are n gradient equations:

∂S ∂β j = 2 ∑ i r i ∂r i ∂β j = 0 (j = 1, ..., n) (2.53)
For a non-linear model, the partial derivatives ∂r i

∂β j are in terms of desired parameters as well as independent variables. Hence, it is usually impossible to get a closed solution. Instead, iterative methods are widely chosen to approximate optimal results in successive manners:

β j ≈ β k+1 j = β k j + ∆β j
(2.54) where k denotes the kth round of iteration. ∆β is known as the shift vector. In each step of iteration, the model f (x, β) can be linearly approximated by a rst-order Taylor series expansion at β k :

f (x i , β) ≈ f (x i , β k ) + ∑ j ∂f (x i , β k ) ∂β j (β j -β k j ) ≈ f (x i , β k ) + ∑ j J ij ∆β j (2.55)
where J ij is the (i, j)th element of Jacobian matrix # J.

J =    ∂f (x 1 , β k )/∂β 1 • • • ∂f (x 1 , β k )/∂β n . . . . . . . . . ∂f (x m , β k )/∂β 1 • • • ∂f (x m , β k )/β n    (2.56) Therefore, J ij = -∂r i
∂β j . and the residuals in Eq. 2.52 are represented as:

r i =y i -(f (x i , β k ) + ∑ j J ij ∆β j ) =∆y i - n ∑ j=1 J ij ∆β j (2.57)
# https://en.wikipedia.org/wiki/Jacobian_matrix Substituting Eq. 2.57 into Eq. 2.53, and rearranging the equations in terms of matrix, we have:

(J T J)∆β = J T ∆y (2.58)
This equation forms the basis for many other iterative approaches of non-linear least squares problem.

Gauss-Newton Algorithm

Gauss-Newton algorithm (GNA) [Bjorck 1996] is very popular to solve non-linear least squares problems. Its biggest advantage is to avoid computing second derivatives.

According to Eq. 2.58, the shift vector for each iteration step can be solved by linear least squares as in Eq. 2.32 ∆β = (J T J) -1 J T ∆y (2.59) Then, in each iteration step, estimated parameters are updated as:

β k+1 = β k + (J T J) -1 J T r(β k ) (2.60)
where the residual r(β k ) is expressed as:

r(β k ) = y -f (x, β k ) (2.61)
In summary, given an initial guess of β 0 at the beginning, Gauss-Newton algorithm iteratively updates the required parameters as in Eq. 2.60 until r(β k ) closes to zero or ∆β k is tiny.

2.3.2.2

Levenberg-Marquardt Algorithm

Levenberg-Marquardt algorithm (LMA) [Levenberg 1944,Marquardt 1963] provides a more robust solution of non-linear least squares than Gauss-Newton algorithm (GNA). In many cases, it converges to local minimum, even if the initial guess starts very far from the nal minimum. Like Gauss-Newton algorithm, Levenberg-Marquardt algorithm is also an iterative method. Recall to Eq. 2.58, which is the basis of Gauss-Newton algorithm. Levenberg [Levenberg 1944] rstly replace this equation by a "damped version":

(J T J + λI)∆β = J T ∆y (2.62)
where I is the identity matrix. λ is a non-negative damping factor, which is adjusted at each iteration. If the S in Eq. 2.51 reduces rapidly, a small value of λ can be used, making the algorithm closer to Gauss-Newton algorithm. Whereas, if an iteration causes insucient reduction of the residual, λ can be increased, giving a step closer to the gradient descent direction (note that the gradient of S with respect to β equals -2(J T ∆y)). Levenberg's algorithm has the disadvantage that if the value of the damping factor λ is large, inverting J T J + λI is not used at all. Marquardt proposed in [Marquardt 1963] to substitute Eq. 2.62 as:

(J T J + λdiag(J T J))∆β = J T ∆y (2.63) where diag(J T J) is a diagonal matrix consisting of the diagonal elements of J T J.

This equation is the basis of Levenberg-Marquardt algorithm. Similar to Gauss-Newton method, after giving an initial guess β 0 , LMA calculates the shift vector as:

∆β = (J T J + λdiag(J T J)) -1 J T ∆y (2.64)
Then, in an iteration step k + 1, parameters are updated as:

β k+1 = β k + (J T J + λdiag(J T J)) -1 J T r(β k ) (2.65)
The iterative procedure stops either r(β k ) is close to zero or ∆β k is tiny.

Robust Estimation methods

Robust estimation methods, or robust statistical methods, is designed to deal with outliers in data. Another motivation is to provide methods with good performance when there is a mixture of models. From statistical theory, M-estimator or least-median squares [Godambe 1991] are proposed. While here, we only describe RANSAC algorithm, a widely used robust estimation method developed from computer vision community. RANSAC (RANdom SAmple Consensus) [START_REF] Fischler | [END_REF]] is a general parameter estimation approach designed to cope with a large proportion of outliers in input data.

It is also a non-deterministic algorithm, which means that it gives a reasonable result by a certain probability.

In aforementioned least squares methods, all the observational data are assumed to be as inliers data t the model only subject to noise at some extent. However, outliers, data do not t the model at all, are pervasive in real applications. The outliers maybe come from extreme values of noise, erroneous measurements, incorrect hypotheses about the interpretation of data, etc.

A simple example is line tting in two dimensions from a set of observations contains both inliers and outliers. As drawn in Fig. 2.13, observational data of a 2D line are severely contaminated by outliers (marked in red in Fig. 2.13 (b)). If we directly estimate the 2D parameters by ordinary least squares, the result would be far from the true model, as shown in Fig. 2.13 (a). While for RANSAC, the outliers have no inuences on the nal results, as drawn in Fig. 2.13 (b). The procedure of RANSAC algorithm is show in Algorithm 1 A trivial problem for RANSAC is how to judge the estimated model. The criterion of the absolute errors of the inliers might make the algorithm pick the model tting a minimal set of inliers (as the model will result in zero error estimation). Usually, the "best model" is typically chosen as the one with the largest consensus set. Therefore, the model with too few satised points is rejected in each step.

The number of iterations N is chosen high enough to ensure that the probability p ransac (usually set to 0.99) that at least one of the sets of random samples contains only inliers. Thus, p ransac gives the probability that the algorithm produces a useful result. Let p inlier represent the probability that any selected data point is an inlier, that is,

p inlier =
number of inliers in data number of points in data (2.66) and p outlier = 1p inlier the probability of observing an outlier. K iterations of minimum L number of points are required, where

1 -p ransac = (1 -p L inlier ) K (2.67)
and taking the logarithm of both sides, leads to

K = log(1 -p ransac ) log(1 -(1 -p outlier ) L ) (2.68) 2.4
Image Local Feature Detectors and Descriptors

Local feature detectors and descriptors are of the most important research areas in computer vision society. Numerous image features have been proposed during the past decades. Detecting specic visual features and associating them across dierent images are applied in a diverse computer vision tasks, including stereo vision [START_REF] Tuytelaars | Wide basebase stereo matching based on local afane invariant regions[END_REF],Matas 2002], vision based SLAM [START_REF] Se | Mobile robot localization and mapping with uncertainty using Scale-Invariant visual landmarks[END_REF]], image stitching [START_REF] Brown | [END_REF][START_REF] Li | [END_REF]] and object recognition [Schmid 1997[START_REF] Sivic | [END_REF]]. According to the types of image feature's usages, there are three main categories:

• Detecting image features for semantic representation. For example, edges detected in aerial images often correspond to roads.

• Stable localizing points through dierent images. What the features actually represent is not really important, as long as their location is accurate and stable enough against changing environments. This maybe the most common usage applied broadly in camera calibration, 3D reconstruction, pose estimation, image registration and mosaicing. A classic example is using Kanade-Lucas-Tomasi(KLT) tracker [Shi 1994] to track image features.

• Image features and their descriptors can be used to interpreted images, such as object recognition, scene classication. In fact, what the feature descriptors really represent are unimportant, the goal is to analyze their statistical performances in an image database.

In this section, we only describe conventional feature detectors and descriptors and their usages in the second category. In essence, localizing a same feature between dierent images is to nd its counterparts. A typical procedure is given as follows:

1. Local invariant feature detector is performed to identify a set of image locations (point, edge or region), which are stable against the variation of environment.

2. A vector carrying on certain visual information around each detected image feature are computed as feature descriptor.

3. Associating detected features between dierent images acquired under various situations. This process could be completed by feature tracking or matching feature descriptors.

2.4.1

Local Invariant Feature Detectors [START_REF] Tuytelaars | [END_REF]] gives a comprehensive survey about local invariant feature detectors. Here, we brief this survey and add some latest progresses on image feature detectors.

A local invariant feature is an image pattern which diers from its neighborhood. It can be in form of a point, an edge or a small image patch. The "invariant" requires relative constant detecting results under certain transformations, such as rotation, scale, illumination changes etc. In general, the following properties are used to evaluate the quality of feature detectors:

• Repeatability: Given two images of the same object or scene acquired under dierent conditions, the percentage of features that occur in both images are dened as repeatability. The repeatability maybe the most important property, which represent the invariance and robustness of a feature detector to various transformations.

• Quantity: The number of detected features should be suciently large, such that a reasonable number of features are detected even on small objects. However, the optimal number of features depends on the application. Ideally, the number of detected features should be adaptable over a large range by a simple and intuitive threshold. The density of features should reect the information content of the image to provide a compact image representation.

• Accuracy: The detected features should be accurately localized, both in image location, as with respect to scale and possibly shape. This property is particularly important in stereo matching, image registration and pose estimation, etc.

• Eciency: Preferably, the detection of features in a new image should allow for time-critical applications.

Since feature detectors are the very fundament of computer vision, numerous approaches have been proposed. In this section, we only introduce several most representative methods.

Corner Detectors

A corner can be dened as a point for which there are at least two dominant and dierent edge directions in its neighborhood. An example is drawn in Fig. 2.14. There are three points (pixels) point a on the corner, point b on the edge and point c inside the object. For the point c, the surrounding pixels (within solid square) remain invariant in all direction. As for the point c, the neighboring pixels change in one direction (perpendicular to the edge). Whereas, the pixels around the point a dier from each other in all directions.

Moravec Corner Detector: It was Moravec [Moravec 1979] who rstly proposed a corner detector nding points that have local maximum in a directional variance measure. The proposed method analyzes a shifting window around a pixel in all directions. As mentioned before, if the pixel is on the corner, none of nearby patches will look similar. The corner strength is dened as the smallest sum of squared differences (SSD) between two patches. If this value is a maximum in local, then, a corner point is found. However, this detector is anisotropic as the intensity variation is only calculated at a discrete set of shifts. Therefore, the operator is not rotationally invariant, which will cause the detector to have poor repeatability rate. Harris Corner Detector: Harris and Stephens [START_REF] Harris | [END_REF]], improved Moravec's method by considering the dierential of the corner score with respect to direction directly, instead of using shifted patches. This corner score is called the autocorrelation matrix, which describes the gradient distribution in a neighborhood of a point. For every pixel (x, y) in an image I, this matrix is expressed as:

M (x, y) = [ ∑ u,v w u.v • [I x (x r , y r )] 2 ∑ u,v w u,v • I x (x r , y r )I y (x r , y r ) ∑ u,v w u,v • I x (x r , y r )I y (x r , y r ) ∑ u,v w u.v • [I y (x r , y r )] 2 ]
(2.69) where I x and I y represent the derivatives of the image in x and y directions respectively, (x r , y r ) = (x + u, y + v), and w(u, v) is a window representing a weighting function, e.g. a binary rectangular window. [START_REF] Harris | [END_REF]] proposed to use a Gaussian window w(u, v) = exp-(u 2 + v 2 )/2σ 2 since it is isotropic. By analyzing the eigenvalues λ 1 and λ 2 of M (x, y), the point can be classied as :

• If λ 1 ≈ 0 and λ 1 ≈ 0, then, this pixel (x, y) has no interest.

• If λ 1 ≈ 0 and λ 2 has large positive value, then, an edge is found.

• If λ 1 and λ 2 both have large positive values, then, a corner is found.

To reduce computational expense, [START_REF] Harris | [END_REF] proposed the following corner score, which depends on the eigenvalues, but avoid direct calculation:

c(x, y) = λ 1 λ 2 -k • (λ 1 + λ 2 ) 2 = det(M (x, y)) -k • trace 2 (M (x, y)) (2.70)
where det( * ) and trace( * ) denote the determinant and trace of a matrix respectively. k is determined empirically (usually in the range of 0.04-0.15). Therefore, Harris corner detector does not actually compute the eigenvalues of Harris matrix in Eq. 2.69. Instead, it calculates the Harris corner score. Shi-Tomasi's "Good Feature to Track": Shi and Tomasi [Shi 1994] theoretically analyzed which feature is good enough to track. Based on the assumption of an ane image transformation, they found that it is more convenient and accurate to use the smallest eigenvalue of the autocorrelation matrix as the corner score.

c(x, y) = min(λ 1 , λ 2 ) > θ • max x,y c(x, y) (2.71)
where θ is a predened percentage factor to control the minimum score. Compared to the Harris score (Eq. 2.70), this requires an additional square root computation on each pixel. Non-Maximum Suppression: Performing corner detectors or other detectors (e.g. object detectors) in a manner of scanning window usually results in multiple responses within high interest regions. This eect makes detecting results too concentrated in a region. A standard method to deal with this is to remove detector responses in the neighborhood of detections with locally maximal score, which is called Non-Maximum Suppression (NMS). For each corner candidate within a region, NMS sets the corner score for this point to zero if its corner score is not larger than the corner score of all candidates within this region. After non-maximum suppression, the corners are the locally maximum valued pixels. An example is shown in Fig. 2.15 Features from Accelerated Segment Test (FAST): Rosten and Drummond [START_REF] Rostern | [END_REF], Rosten 2006] developed a high-speed corner detector: Features from Accelerated Segment Test, coined as FAST. This method performs on a discretized circle (Bresenham circle) around a candidate point p as shown in Fig. 2.16. The black point in the center is the candidate point p. The 16 grey points are the discretized approximation of the outlined circle around it.

The FAST detector classies p as a corner if there is a large number of pixels on a circle of xed radius around p such that these pixels are all signicantly brighter, or darker, than p by a threshold t. The algorithm is further accelerated by a trained decision tree to operate as few pixels as possible. With this decision tree, only 2.26 pixels in average are tested [Rosten 2006]. To perform a Non-Maximum Suppression, Figure 2.16: Bresenham circle used in FAST detector [Rosten 2006] the following corner score is computed for each candidate point:

c(p) = max{ ∑ q∈S + |I q -I p | -t, ∑ q∈S - |I q -I p | -t} (2.72)
where S + is the subset of pixels on the circle that are brighter than p (by the threshold t) and S -is the subset of pixels that are darker than p (by the threshold t).

Interest Point Detectors

Interest point detector is an extension of corner point. A point can be characterized as an interest point if:

• The local visual information (e.g. texture, gradient) around the point is rich.

• Its position is stable under certain image transformations, such as scale/illumination/ane changes.

A corner point detector (described in last section) is a special kind of interest point detectors that it is mainly designed for robust image feature tracking. Here, we discuss more general cases of detecting interest points (or regions) for more general tasks. Please note that in the history of computer vision, the applications of terminologies "blob detector" and "interest point detector" are overlapping. Output of blob detectors are usually well-dened point positions, which may correspond to local extremas in an operation applied in local regions. (2.73) with L xx , L xy and L yy refer to the second order derivatives of a Gaussian smoothed image at (x, y). Gaussian smoothing is performed due to second order derivative's extreme sensibility to noise:

L(x, y; σ) = G(x, y, σ) * I(x, y) (2.74)
where G(x, y.σ) is a Gaussian kernel at scale σ:

G(x, y, σ) = 1 2πσ 2 e -(x 2 +y 2 )/(2σ 2 )
(2.75) Hessian matrix encodes rich information about local image structure.

Laplacian of Gaussians (LoG) is probably the rst commonly used blob feature detector. It is dened by the trace of Hessian matrix: ∇ 2 L = L xx + L yy . LoG preserves the edges and corners on an image, as seen in Fig. 2.17. These edges and corners are good candidates for a good feature detector. However, the LoG operator strongly depends on the size of blob in the image and the size of Gaussian kernel used for smoothing. In order to automatically detect blobs of dierent scale in the image, a multi-scale approach is proposed. A straightforward approach is to use the scale-normalized Laplacian operator (normalized LoG):

∇ 2 norm L(x, y; σ 2 ) = σ 2 (L xx + L yy ).
(2.76) This operator is to detect scale-space extremas, which are local maxima/minima of ∇ 2 norm L(x, y; σ 2 ) with respect to both space and scale changes [Lindeberg 1994[START_REF] Bretzner | [END_REF]]. In practice, computing second order derivatives is computationally intensive. Several approximations described in the following (Dierence of Gaussians and Fast Hessian) are proposed to speed up the computation time.

Dierence of Gaussians (DoG): The Dierence of Gaussians detector (DoG), as a part of SIFT [Lowe 2004], is proposed to approximate normalized LoG. Specically, a DoG operator is given by

D(x, y, σ) = L(x, y, λσ) -L(x, y, λσ) = (G(x, y, kσ) -G(x, y, λσ)) * I(x, y) (2.77)
where G(x, y, σ) is the Gaussian kernel dened in Eq. 2.75. λ is a constant multiplicative factor. The convolution results L(x, y, kσ) with dierent k create a scale space. SIFT uses octaves of scale spaces to nding extremas of DoG. The approach to quickly construct D(x, y, σ) in several scale spaces is shown in Fig. 2.18.

For each octave of scale space, the initial image is repeately convolved with Gaussians to produce the set of scale space images shown on the left of Fig. 2.18. Adjacent Gaussian images are subtracted to produce the dierence-of-Gaussian images on the right of Fig. 2.18. After processing each octave, the Gaussian image is down-sampled by a factor of 2, and the process repeated. Key points are then extracted at local minima/maxima of the DoG images through scales.

The DoG approximation greatly reduces the computation cost it can be computed by a simple image subtraction. As an approximation of normalized LoG operator, DoG is scale invariant.

Fast Hessian Detector Fast Hessian detector, proposed by [START_REF] Bay | [END_REF]] as a basis for SURF (Speeded Up Robust Features), is based on the determinant of Hessian matrix with scale factor σ:

det H(x, y, σ) = L xx (x, y, σ) • L yy (x, y, σ) -L 2 xy (x, y, σ) (2.78)
To further accelerate the computation of Gaussian second order derivatives, [START_REF] Bay | [END_REF]] approximates them with simple box lters (Fig.

det[H(x, y, σ)] ≈ D xx (σ) • D yy (σ) -(0.9D xy (σ)) 2 (2.79)
Eq. 2.79 is taken as a measure of key point, and performed in the image. The responses are stored in a response map, where local extremas are found and rened.

In [START_REF] Bay | [END_REF]], it reports that SURF is more than ve times faster than DoG.

Center-Surround Extrema (CenSurE): CenSurE [START_REF] Agrawal | [END_REF]] is proposed to execute the task of matching two images in real time for camera motion estimation, especially in dicult environment where there is large image motion between frames.

While SIFT [Lowe 2004] approximates LoG with dierence of Gaussians, CenSurE adopts a further approximation, utilizing bi-level center-surround lters, i. The computation of bi-level lters is even further speeded up by slanted integral image, as a modication of original integral image [Viola 2001]. Following the output of a bi-level lter, a 3 × 3 × 3 neighborhood non-maximum suppression and an edge/line response suppression are performed to rene the results. The authors claimed that CenSurE outperforms the other detectors (e.g. SIFT, SURF, FAST, Harris corner) in stability and accuracy.

Feature descriptors

After extracting interest points, a descriptor is created to identify and match them between dierent images. The most straightforward description for a feature point is the intensity appearance within a patch around the point itself. Similarity metrics, e.g. the sum of absolute dierences (SAD), the sum of squared dierences (SSDs), or the normalized cross correlation (NCC) can be used:

M SAD = ∑ (i,j)∈p |I 1 (i, j) -I 2 (x + i, y + j)| M SSD = ∑ (i,j)∈p (I 1 (i, j) -I 2 (x + i, y + j)) 2 M N CC = ∑ (i,j∈p) I 1 (i, j)I 2 (x + i, y + j) √ ∑ (i,j)∈p I 2 1 (i, j) ∑ (i,j)∈p I 2 2 (x + i, yj) (2.80)
In many cases, local appearance around a feature point is not enough, especially when encountering large changes in orientation, scale or viewpoint. In fact, none of the SSD, SAD and NCC are invariant to any of these changes. Therefore, their use is limited to images taken at nearby positions.

Exploring feature descriptors stable in various changes has been conducted for decades in the computer vision eld. In the following, we will introduce several popular feature descriptors. SIFT descriptor: One of the most popular descriptor for point features is the SIFT [Lowe 2004]. Fig. 2.21 illustrates how the descriptor is computed. At rst, SIFT descriptor is proved to be stable against changes in illumination, rotation, and scale, and even up to 60 changes in viewpoint. Although proposed in 2004, it still among the top-performanced feature points [Gauglitz 2011]. In general, the SIFT descriptor can be computed for corner or blob detectors. However, its performance will decrease on corners because that corners are not as distinctive as for blob detectors.

SURF descriptor: Although SIFT achieve great success, its 128-dimensional feature descriptor costs intensive computation in feature matching or indexing. SURF descriptor [START_REF] Bay | [END_REF]] simplies SIFT descriptor while keeping almost the same matching performance. [START_REF] Bay | [END_REF]] rstly computes an orientation for each keypoint. Haar wavelet lters (see in the left gure of Fig. 2.22 (a)) in x and y directions are convoluted with a circular region around the keypoint. Similar to SIFT, the lter outputs are weighted by a Gaussian and then represented as 2-dimensional vectors, which are summed up with a rotating angular window. The longest resulting vector determines the orientation of the keypoint. Then, a square region centered at the keypoint, and

v = [ ∑ d x , ∑ d y , ∑ |d x |, ∑ |d y |]
(2.81)

where d x and d y are responses of Haar wavelet lter in x and y directions. Hence, this results in a descriptor vector for all 4 × 4 subregions of a length equals to 64. This process is illustrated in Fig. 2.22 (b).

ORB descriptor: In recent years, several binary descriptors have been designed to achieve high computation speed to meet requirements of mobile visual computing. Here, we will introduce two typical binary descriptors: Oriented FAST and Rotated BRIEF (ORB) [Rublee 2011] and Binary Robust Invariant Scalable Keypoints (BRISK) [Leutenegger 2011].

ORB feature [Rublee 2011] can be viewed as an improved combination of FAST and BRIEF [Calonder 2010] such that it is based on an oriented FAST corner detector plus a rotated BRIEF descriptor. The detector component relies on intensity centroid [Rosin 1999] to compute orientation for every detected FAST feature point. As for the feature descriptor, the BRIEF descriptor is a bit string description of an image patch constructed from a set of binary intensity tests. Consider a smoothed image path, p. A binary test τ is dened by:

τ (x, y; p) = { 1 : p(x) < p(y) 0 : p(x) > p(y)
where p(x) is the intensity at point x. The feature is dened as a vector of n binary tests:

f n (p) = ∑ 1≤i≤n 2 i-1 τ (p; x i , y i ) (2.82)
In order to let BRIEF be invariant to in-plane rotation, [Rublee 2011] proposed to steer BRIEF by the orientation of keypoints. For any feature set of n binary tests at location (x i , y i ), dene the 2 × n matrix:

S = ( x 1 • • • x n y 1 • • • y n ) (2.83)
Using the patch orientation θ estimated by intensity centroid, and the corresponding rotation matrix R θ , we get a "steered" version S θ = R θ S. The steered BRIEF operator becomes:

g n (p, θ) := f n (p)|(x i , y i ) ∈ S θ (2.84)
To recover from the loss of variance in steered BRIEF and to reduce correlation among the binary tests, [Rublee 2011] develops a learning method for choosing a good subset of binary tests. The result is called rBRIEF, which signicantly improves the variance and correlation of steered BRIEF. The rBRIEF descriptor is taken as ORB's descriptor. The authors demonstrate that ORB is at two orders of magnitude faster than SIFT, while performing as well in many situations.

BRISK descriptor: Another binary descriptor is BRISK, which is proposed in [Leutenegger 2011]. The detector of BRISK is an extension of FAST detector by searching maxima not only in the image plane, but also in scale-space using the FAST score in Eq. 2.72 as a measure for salience. Similar to ORB, BRISK descriptor is also composed as a binary string by concatenating the results of simple brightness comparison tests. The direction of each keypoint is also computed to achieve rotation invariance. Furthermore, BRISK descriptor selects the intensity comparisons maximizing descriptiveness between different points. The key concept of BRISK descriptor is using a pattern (as illustrated Figure 2.23: BRISK sample pattern [Leutenegger 2011] in Fig. 2.23) to sample the neighborhood of the keypoint. The pattern denes N sampling locations equally spaced on concentric circles with the keypoint. Let us consider one of the N •(N -1)/2 sampling-point pairs (p i , p j ). The smoothed intensity values at these points which are I(p i , σ i ) and I(p j , σ j ) (The Gaussian smoother in BRISK is dierent to the position of p). The local gradient g(p i , p j ) is:

g(p i , p j ) = (p j -p i ) • I(p j , σ j ) -I(p i , σ i ) ∥p j -p i ∥ 2 (2.85)
Considering the set A of all sampling-point pairs, the authors dene a subset S of short-distance pairing and another subset L of long-distance pairings:

S = {(p i , p j ) ∈ A| ||p j -p i || < σ max } ⊆ A L = {(p i , p j ) ∈ A| ||p j -p i || > σ min } ⊆ A (2.86)
The pattern direction is estimated through the point pairs in L:

g = ( g x g y ) = 1 L • ∑ p i ,p j ∈L g(p i , p j ) (2.87)
where L is the number of point pairs. In order to make a rotation and scale normalized descriptor, BRISK rotates the sampling pattern by α = arctan 2(g x .g y ) around a keypoint. The bit-string descriptor is created by comparing all the short distance intensity comparisons:

b = { 1 I(p α j , σ j ) > I(p α j , σ i ) (p α i , p α j ) ∈ S 0 otherwise
The experiments in [Leutenegger 2011] reveal BRISK's high quality performance as in state-of-art algorithms, with a lower computational cost.

Associating feature points through images

As described in the beginning of Sec. 2.4, one of the most important need when using image feature points is to stable localizing the same point through dierent image. Many computer vision applications, e.g. stereo matching, pose estimation, augmented reality and object recognition, are relied on precise feature point associating. In this section, we introduce two approaches to associate feature points: the rst is feature tracking, the second is feature matching. KLT tracker is a classic image feature point tracking algorithm in computer vision area. It is based on the early work of [Lucas 1981], then is developed fully by [START_REF] Tomasi | [END_REF]], and is explained in details in [Shi 1994]. Later, [Bouguet 2000] gives a pyramidal implementation to cope with large displacement. Let two images I(x) = I(x, y), J(x) = J(x, y) refer to the rst image and second one, respectively. For an image point u = [u x , u y ] on the rst image I, KLT tracker is to nd its corresponding location

v = u + d = [u x + d x , u y + d y ]
T on the second image J. Let W(x; p) be the warp from image I to J, where p is a vector of parameters. The best alignment minimizes image dissimilarity:

∑ [I(W(x; p)) -J(x)] 2
(2.88) Eq. 2.88 expresses nonlinear optimization. Therefore, an iterative method is used to solve p. Assume that a p is known and the best increment ∆p is sought that the problem is modied as minimizing the following quantities:

∑ x [I(W(x; p + ∆p)) -J(x)]
(2.89)

The above problem is iteratively solved with respect to ∆p: p ← p + ∆p. Eq. 2.89 can be linearized by rst order Taylor expansion:

∑ x [I(W(x; p)) + ∇I ∂W ∂p ∆p -J(x)] 2
(2.90)

∇I = [ ∂I ∂x , ∂I
∂y ] is the gradient image computed at W(x; p). The term ∂W ∂p is the Jacobian of the warp. Computing the derivative of Eq. 2.90 with respect to ∆p, we obtain:

2 ∑ x [∇I ∂W ∂p ] T [I(W(x; p)) + ∇I ∂W ∂p ∆p -J(x)]
(2.91)

Setting Eq. 2.91 to zero yields:

∇p = H -1 ∑ x [∇I ∂W ∂p ] T [J(x) -I(W(x; p))] (2.92)
where H is the Hessian matrix:

H = ∑ x [∇I ∂W ∂p ] T [∇I ∂W ∂p ]
(2.93) By putting the above equations in an iterative approach, we get a whole description of KLT tracker, as summerized by Algorithm 2 Warp I with W(x; p)

3:
Warp the gradient ∇I with W(x; p)

4:

Evaluate the Jacobian ∂W ∂p and compute the steepest descent image ∇I ∂W ∂p 5:

Compute the Hessian H according to Eq. 2.93 6:

Compute ∆p according to Eq. 2.92 7:

Update the parameters p ← p + ∆p

8: until ||∆p|| < ε 2.4.3.2 Feature matching
Similar to feature tracking, the process of feature matching is also for nding corresponding features in dierent images. The dierence is that, feature matching is to nd correspondences between existed feature points, feature tracking is to compute a feature point's new location in a new image. Feature matching is achieved by comparing feature descriptors between dierent feature points. Usually, a similarity metric, such as SSD, NCC, or Euclidean distance between two feature vectors, is computed for searching correspondences.

In essence, the feature matching problem is a Nearest Neighbor (NN) search problem, or Post Oce Problem as: Denition 1. Given a set P ∈ R d of n points and a query point q ∈ R d , the nearest neighbor search of q is:

N N (q) = arg min p∈P |p -q| (2.94)
where q is a feature point in d-dimensional feature space. The most simple solution of NN search is to compare all the candidate points with the query point q, nd a point with the highest similarity score or the closest distance.

In feature matching, simple NN search may result with features in the second image matching with more than one feature in the rst image. A mutual consistency check is usually used to eliminate this eect by reversely comparing the features in the second image with features in the rst image. Only pairs of corresponding features matched in both directions are accepted as correct matching. However, disadvantage of exhaustive NN search is obvious: its computation grows in a quadrical speed when the number of candidates increases. To speed up the searching process, several indexing structures, such as multidimensional search tree (e.g. k-d tree [START_REF] Friedman | [END_REF]) or hash table [START_REF] Lamdan | [END_REF]], are proposed to rapidly search for features near a given feature. An example of 3-dimensional k-d tree is shown in Fig. 2.24. It can achieve eciently NN searching by partitioning feature spaces into indexed leaves. An another feature matching strategy is to search for potential correspondences in an expecting region in the second image. The size of regions can be predened or predicted by the motion model of feature points. More works can be refered in [Chili 2008[START_REF] Handa | [END_REF]].

Conclusion

In this chapter, we rstly describe the stereo vision measurement model. Since our stereo vision system consists of two quasi-parallel monocular camera, the stereo vision measurement model is based on monocular camera pinhole model. The distortion model of the lens is also considered. The stereo vision model is very important and forms the basis of many further application in the following chapters.

In the second, we review second linear/non-linear least squares methods. These methods are commonly used to solve estimation problems. For linear least square methods, ordinary and total least squares methods are introduced. Ordinary least squares are widely used in many computer vision problems, while the total least squares is an improved version of least squares when considering all the errors. For non-linear least squares, Gaussian-Newton method and Levenberg-Marquardt algorithm are described. We will use these two methods to estimate parameters in the following chapters.

At last, we present an overview of feature detectors/descriptors. In literature, a lot of feature detectors/descriptors are proposed because of wide applications. The number of proposed feature detectors/descriptors is so big that many researchers often feel confused. We try to summarize the development of feature detectors/descriptors by dierent categories. In the following chapter, we will analyze their performances in visual odometry. In this chapter, we provide a stereo vision based method to estimate the egomotion of a moving intelligent vehicle in urban environments. Estimation of selfmotion is important for an intelligent vehicle since:

• It provides a complementary localizing approach for GPS. In urban environments, GPS signals could disappear or be disturbed due to buildings, tunnels, etc. Without any prior knowledge of the environment nor a predened motion model of the vehicle, vision based ego-motion estimation can estimate the path of a camera-equipped vehicle by calculating the ego-motion between consecutive images in a video ow. Hence, vision based ego-motion estimation could be used for vehicle localization when GPS signals are not available.

• It is the fundament of further analyzing of a dynamic environment. When the intelligent vehicle is driving in a dynamic environment, understanding selfmotion is the rst step before performing independent moving object detection, tracking and recognition.

In this chapter, we present a stereo vision based method to estimate the ego-motion of a moving intelligent vehicle. The contributions are comparing tracking based approaches with matching based approaches, as well as performances dierent feature detectors. The objective is to select the best feature point detector and the best circular point association for urban environment perception.

Introduction

Vision based ego-motion estimation, or, visual odometry (VO), computes movement of a vehicle by the outputs of on-board single or multiple cameras. Usually, it is used in the eld of robotics or intelligent vehicles. Visual odometry operates by incrementally estimating the pose of the vehicle through the changes of images caused by the motion. Compared to global positioning system (GPS), visual odometry could work in more general environments, such as indoor, underwater and even in the planet Mars [START_REF] Matthies | [END_REF]]. Advantages of visual odometry with respect to the other navigation sensors, e.g. inertial measurement units (IMU), laser range nder, is that it is not aected by wheel slip and it provides more accurate trajectory estimates in a low expense.

In literature, visual odometry (VO) is a particular case of structure from motion (SFM) [Dellaert 2000]. For an image sequence captured by a moving camera, SFM not only recovers the camera's trajectory, but also reconstructs the 3D structures of all the vision measurements. While in many applications, full 3D reconstruction is not necessary. Visual odometry focuses only on accurately estimating 3D motion of the camera frame by frame.

Approaches to visual odometry range from dense optical ow [START_REF] Mccarthy | [END_REF]], matching sparse salient image regions in consecutive images, to semiparametric regression from sparse optical ow [Guizilini 2012] proposed recently. In most navigation/localization applications, real-time performance is a requirement. Hence, sparse correspondence based methods are usually preferable to the other kinds of methods. It was Moravec [Moravec 1980] who rstly proposed a vision-based method to estimate a planetary rover's ego-motion. The proposed method is based on tracking sparse salient points detected in stereo image pairs. Similar early researches ( [START_REF] Matthies | [END_REF],Matthies 1989[START_REF] Olson | [END_REF]) in visual odometry were driven by NASA Mars exploration. [Moravec 1980]'s method rstly detects corner points (Moravec corner point in Sec. 2.4.1.1) in the images of a sliding camera. Then, 3D positions of the corner points are computed by triangulation. Finally, rigid motion information is estimated by aligning the triangulated 3D points between two consecutive frames. As a summary, [Moravec 1980] proposed a pipeline to solve visual odometry problem based on sparse feature points its main functioning parts remain unchanged today. The framework of a typical visual odometry is drawn in Fig. 3.1. At rst, certain kind of salient feature points are detected and then associated through consecutive input images. Based on the associated pairs of feature points, 2D-2D/3D-2D/3D-3D constraints can be utilized to estimate moving information. Sometimes, a local optimization method (e.g. bundle adjustment [START_REF] Bill | [END_REF]) can be optionally chosen to improve the results.

Under the framework in Fig. 3.1, many researches have been conducted to improve visual odometry's accuracy and robustness. Since a binocular stereoscopic system is mounted in our platform, we only review the stereo vision based visual odometry methods. [START_REF] Matthies | [END_REF],Matthies 1989] tain more accurate trajectory results, when compared with those in [Moravec 1980]. [START_REF] Lacroix | [END_REF]] detects candidate key points from dense stereo by analyzing the correlation function around its peaks. This kind of key point demonstrates advantages in accurate feature localization. Later works, as [START_REF] Cheng | [END_REF], Cheng 2007], utilize a similar curvature of the correlation around Harris corner point to dene an error covariance matrix. Furthermore, they adopt RANSAC in least squares estimation for rejecting outliers. For each stereo pair, they reconstruct the 3D positions of feature points by triangulation. The motion estimation is treated as a 3D-3D point alignment problem. To build correspondences of feature points through frames, [Nister 2004] uses feature matching to replace previous feature tracking. It could ameliorate feature drift problem when tracking features through multiple frames. The motion information is computed from 3D-2D point constraints. A similar procedure is proposed in [START_REF] Geiger | [END_REF]]. It greatly reduces the feature matching time by extracting binary features. Another dierent approach proposed in [START_REF] Comport | [END_REF]] estimates the motion from 2D-2D image point matches without their 3D positions. This solution derives from a "quadrifocal tensor". The authors claim that more accurate motion estimation could be reached by the 2D-2D based motion estimation.

In the following, we will follow the framework of classic stereo visual odometry to estimate vehicle's ego-motion. The contribution of our works is the comparison of various feature detectors and the dierent approaches of circular feature association.

Circular Feature Points Detection and Association

The objectives of this section are nding local invariant salient points in successive stereo image pairs and establishing their mutual correspondences. As illustrated in Fig. 3.1, these steps provide necessary information for the following motion estimation. In our practice, we detect and associate the feature points in a circular manner as in [START_REF] Geiger | [END_REF]]. For an acquired stereo image pair at time t -1 and t, circular point-to-point correspondences are established in all the four images: current left/right images and previous left/right images. This is depicted in Fig. 3.2. In this step, a feature point (u t l , v t l ) detected in the current left image at time t is associated with corresponding locations in the other images as:

(u t r , v t r ), (u t-1 l , v t-1 l ), (u t-1 r , v t-1 r ).
Therefore, the four positions of a same point in four images form a circle: In Sec. 2.4, we have reviewed several commonly used image local detectors, descriptors, as well as tracking and matching approaches. To seek a circular pointto-point correspondences between successive images, two distinct approaches are possible:

(u t l , v t l ) ↔ (u t r , v t r ) ↔ (u t-1 l , v t-1 l ) ↔ (u t-1 r , v t-1 r ) t-
• The rst approach is to detect feature points in one image, then, track the detected points in other images.

• The second approach is to detect feature points in all the four images, then, match these points by their feature descriptors. is within a tolerance range, the circular tracking is succeeded. Otherwise, the KLT tracker failed to track this point within a circle. An example of KLT tracker based circular point association is demonstrated in Fig. 3.3 (a). The yellow circle shows a successfully circularly tracked feature point p t l,1 from the current left image to the other three images. The green circle shows an unsuccessfully tracked circle due to unacceptable displacement.

Circular Feature Point Association by Matching

As for feature matching based circular feature point association, feature points and their descriptors are computed separately in all the four images. Given a feature point p t l in the current left image, we search its correspondences in the other three images by comparing their descriptors. To avoid exhaustive search in the whole image, we dene a region of interest (ROI) to speed-up the matching process. Since the points are matched between consecutive frames, the change between images is usually not too large. Hence, the ROI is set as a xed rectangular area around p t l in all the other three images. Similarly, the matching process is rened by comparing the displacement of matched feature points. The matched points that have large distance between each other are discarded to assure the quality of matching. A typical example is depicted in Fig. 3.3 (b). The yellow circle shows a successfully circularly matched feature point p t l,1 from the current left image to the other three images. The green circle shows an unsuccessfully circularly matched feature point p t l,2 due to deviation. The purple dotted rectangle is the ROI to simplify searching process.

In addition, two real examples of circular feature points association by KLT tracker and feature matching are shown in Fig. 3.4 (a) and (b), respectively. represents a successfully associated image feature point and its correspondences:

s t i = {p t l,i , p t r,i , p t-1 l,i , p t-1
r,i }, and N denotes the total number of such subsets. Referring to Sec. 3.1, given a set S t of circularly associated feature points, there are three types of geometrical constraints (3D-3D/3D-2D/2D-2D) [Scaramuzza 2011] to compute ego-motion information. In this section, we only describe how to use 3D-2D constraint to estimate ego-motion.

3D-2D Constraint based Ego-motion Estimation

Suppose a binocular stereo rig is already calibrated and the stereo images are all well rectied, then, real 3D positions of the sparse feature points can be easily reconstructed by triangulation (Eq. 2.23). Taking the left image frame as the reference coordinate system, a subset of associated feature points s t = {p t l , p t r , p t-1 l , p t-1 r } can reconstruct two 3D positions of a same feature point in consecutive times:

P t =   X t Y t Z t   =   (u t l -c u ) • b/∆ t (v t l -c v ) • b/∆ t f • b/∆ t   , P t-1 =   X t-1 Y t-1 Z t-1   =   (u t-1 l -c u ) • b/∆ t-1 (v t-1 l -c v ) • b/∆ t-1 f • b/∆ t-1   (3.1) where ∆ t = |u t l -u t r | and ∆ t-1 = |u t-1 l -u t-1
r | are the disparities at time t and t-1 respectively. b, f and (c u , c v ) correspond to the baseline length, focal length, the position of principal point, respectively. Suppose the displacement between P t-1 and P t is entirely caused by the self-motion of the stereo rig, which amounts to assume the point is static. Then, P t-1 and P t are related by a rigid 3D transformation:

P t = [ R t t-1 T t t-1 0 1 ] P t-1 = M t t-1 • P t-1 (3.2)
where • denotes homogeneous coordinates, R t t-1 and T t t-1 represent the 3×3 rotation matrix and 3 × 1 translation vector from time t -1 to t, respectively. M is a 4 × 4 augmented transformation matrix. Directly computing motion parameters from Eq. 3.2 is based on 3D-3D constraint. However, issuing from the inherent inaccuracy of 3D triangulation, this approach usually results large deviation. To make a step forward, in ideal situation, considering corresponding image position p t , we can have:

pt = K[I|0]M t t-1 • P t-1 (3.3)
where K is the 3 × 3 intrinsic matrix in Eq. 2.12. Since the errors are inevitable in real problems, given a set S t containing N circularly associated feature points, ego-motion can be solved by minimizing:

arg min R,T N ∑ i=1 ||p t l -K[I|0]M t t-1 P t-1 || 2 (3.4)
The minimization of Eq. 3.4 is a typical non-linear least square problem. We use Gauss-Newton method (Sec. 2.3.2.1) to iteratively minimize Eq. 3.4 with respect to the transformation parameters (R, T).

Ego-motion estimation by Gauss-Newton method: In fact, the desired rotation matrix can be represented by three parameters: the yaw angle α, pitch angle β and roll angle γ:

R(α, β, γ) =   cos β cos γ -cos β sin γ sin β sin α sin β cos γ + cos α sin γ -sin α sin β sin γ + cos α cos γ -sin α cos β -cos α sin β cos γ + sin α sin γ cos α sin β sin γ + sin α cos γ cos α cos β   (3.5)
and the translation vector is:

T = [t x , t y , t z ] T (3.6)
Accordingly, Eq. 3.3 can be expanded by substituting Eq. 3.5 and Eq. 3.6 into Eq. 3.3 as:

[ u t v t ] = [ f • X t Z t f • Y t Z t ] (3.7)
where

X t = cos β cos γ • X t-1 -cos β sin γ • Y t-1 + sin βZ t-1 + t x Y t = (sin α sin β cos γ + cos α sin γ) • X t-1 -(sin α sin β sin γ + cos α cos γ) • Y t-1 -sin α cos β • Z t-1 + t y Z t = (-cos α sin β cos γ + sin α sin γ) • X t-1 + (cos α sin β sin γ + sin α cos γ) • Y t-1 + cos α cos β • Z t-1 + t z (3.8)
Therefore, for one data point, the corresponding Jacobian matrix is: The partial derivatives with respect to α, β, γ, t x , t y , t z can be easily deducted from Eq. 3.9. Given N circularly matched feature points, the nal Jacobian matrix can be created by packing all corresponding Eq. 3.9 together. Then, Gauss-Newton method is performed: Given an initial guess of the rotation matrix and translation vector, the residual is calculated according to Eq. 3.4. Next, a shift vector to update the parameters is computed by Eq. 2.59. When the parameters are updated, the residual according to the new parameters is re-computed. Ego-motion parameters are obtained by repeating such steps until the residual reaches a minimum or the maximum number of iterations is reached. In practice, we note that even if we initialize R and T to 0, a couple of iterations are sucient for convergence.

However, in a real dynamic urban environment, the components of circularly associated feature points have three categories: a portion may belong to static surroundings, several points could be from either independent moving objects (moving pedestrians, vehicles, etc.) or inaccurate associations. For the ego-motion estimation, the last two types of associations are outliers to be removed. To be robust against outliers, we wrap previous estimation process into a RANSAC frame (Sec. 2.3.3.1). The RANSAC based stereo visual odometry is shown in Algorithm 3. Randomly select at least 3 observation subsets from S t (three is the minimum number for solving 6 motion parameters).

3:

Minimize the re-projection error in Eq. 3.4 and get new motion parameters by Gauss-Newton method.

4:

Classify inliers according to the new motion estimation.

5:

Update motion estimation if more inliers are obtained. 6: until The reprojection error converge or the number of iteration reaches the maximum 7: Final optimizing in all the inliers to rene the estimated parameters 3.4 Experiments of Stereo Visual Odometry

The described ego-motion estimation approach is evaluated by our experimental vehicle SetCar introduced in Sec. 1.3. The installed stereo vision system (Bumblebee XB3) observes surroundings by stereo image pairs (with a resolution of 1280 × 960) in a frame-rate of 13 -15f ps with a baseline 0.24 of meters. The whole framework is implemented in C++, based on OpenCV library , without any acceleration technique. A laptop with a CPU Intel i7-3770 quad cores 3.40GHZ running in Linux is used to run the software.

To begin with, we rstly evaluate the performances of dierent approaches for circularly associating feature points. We mainly compare tracking based approaches with matching based approaches, as well as performances of dierent feature detectors. The best feature point detector and the best approach of circular point association are selected for nal real experiments in urban environments.

Comparing Dierent Feature Association Approaches

In Sec. 2.4, we have reviewed several image feature detectors and descriptors.

Here, we evaluate their performances in circular feature association. The involved feature detectors are: Good feature to track (GFTT) [Shi 1994], FAST [START_REF] Rostern | [END_REF], Rosten 2006], SIFT [Lowe 2004], SURF [START_REF] Bay | [END_REF]], ORB [Rublee 2011], CenSurE (STAR in OpenCV) [START_REF] Agrawal | [END_REF]], BRISK [Leutenegger 2011], and feature descriptors are: SIFT [Lowe 2004], SURF [START_REF] Bay | [END_REF]], ORB [Rublee 2011], BRISK [Leutenegger 2011]. For a fair comparison, all the feature detectors and descriptors use tuned parameters implemented in OpenCV. Evaluating feature detectors: In [START_REF] Tuytelaars | [END_REF][START_REF] Mikolajczyk | [END_REF], several general criterions for evaluating feature detectors are given (refer to Sec. 2.4). Different to the application of object recognition, which usually has prominent scale or viewpoint changes, in our application, the detectors are used to establish point-topoint correspondences between consecutive video frames. Since the movement of our platform is stable and the frame rate of our stereo-vision system is around 15f ps, the scale and viewpoint changes are not large in our application. Hence, the performances of dierent feature detectors are evaluated from 3 aspects: repeatability, uniformity and speed :

• Repeatability: Given two images of the same scene acquired in dierent conditions, the percentage of features that could be found in both images is dened as repeatability. Repeatability score is dene as:

S r = f - f * (3.10)
where f -is the number of features found in both images, f * is the number of features detected in one of the two images, noted as reference image.

http://opencv.org

• Uniformity: To better estimate motion in dynamic scenes, the detected features should be distributed as even as possible. To evaluate uniformity, we divide the 1280 * 960 size image into n u disjoint identical grids. If the number of feature points located in a grid is more than a threshold, it would be marked as "lled". The uniformity score is:

S u = n - u /n u (3.11)
where n - u is the number of "lled" grids. • Speed: The processing speed is measured in millisecond (ms). Usually, the detectors with high processing speed is more preferable in real-time tasks, such as ADAS in intelligent vehicle applications.

To evaluate repeatability, we convert an image by a homography transformation, in a similar way as in [START_REF] Tuytelaars | [END_REF]]. Then, feature detectors are independently performed in the two images. Next, we wrap the detected feature points in the rst image to the second one by the known homography matrix. A feature point is repeated, if the wrapped position in the second image also has a detected feature point. In our experiments, randomly generated homography transformations conned in a moderate range are applied to 100 images captured by our stereoscopic system in the city of Belfort. Repeatability scores are computed as in Eq. 3.10. The performance of all the measured detectors are shown in Fig. 3.5. In this eval- uation, STAR performs the best, with a score of 0.9, followed by SIFT (0.82) and GFTT (0.8). FAST performs the worst, with a score of 0.64. The results reveal that the DoG (Dierence of Gaussian) based methods (STAR, SIFT) attain the best repeatability score.

GFTT FAST SIFT SURF ORB STAR BRISK

Repeatability Performance

To evaluate uniformity score, each test image is divided into 192 squares with a size of 80 × 80 pixels and the threshold is set to 2. We tested all the 100 selected images and the evaluation results are shown in Fig. 3.6 (a). It can be seen that FAST performs the best with a uniformity score of 0.67, followed by STAR, SURF, GFTT, SIFT and ORB performs the worst. The main reason that FAST reaches the best score is of its huge number of detected feature points. The results of processing FAST is no doubt much more faster than the other detectors thanks to its extremely simplied computation. GFTT, ORB and STAR are almost 10 times slower than FAST. SIFT is the slowest due to the computation of DoG and searching maxima in multi-scale spaces.

Evaluating circular feature point association: As we stated in Sec. 3.2.2, there are two approaches, the KLT tracker based and feature matching based, for circularly associating feature points in four images.

For the matching based approach, feature descriptors can be combined with almost all kinds of detectors, such as "GFTT detector + SIFT descriptor". To avoid tedious, reduplicative and similar comparison, we only take original detectordescriptor combinations into consideration, (e.g. "SIFT detector + SIFT descriptor", "SURF detector + SURF descriptor", "ORB detector + ORB descriptor", "BRISK detector + BRISK descriptor" ). For the KLT tracker based approach, we evaluate following detectors: GFTT, FAST and STAR. Also, we propose two criterions to evaluate their performances:

• Survival rate : Both KLT tracker and matching based methods would eliminate a part of detected feature points due to inaccurate tracking or matching. When tracking or matching feature points through four images, such phenomenon is amplied. The "survived" associated feature points maybe too few to be utilized. Hence, we use survival rate S s to describe the performance:

S s = n lef t n total (3.12)
where n total is the number of all detected feature points in a reference image and n lef t is the number of successfully tracked or matched feature points in the other images.

• Accuracy: As we have described in Algorithm 3, a RANSAC scheme is used to exclude outliers caused by independent moving objects and imprecisely associated feature points. Although RANSAC is excellent enough to cope with outliers, we still hope that, when moving in a static environment the portion of outliers is as small as possible. Therefore, we introduce an accuracy factor: the ratio of the inliers to number the whole feature points.

S a = n inlier n total (3.13)
In our evaluation, the images without independent moving objects are used.

Except the proposed two criterions, the performance of speed is also very important.

• Speed: For KLT tracker based method, the processing time is counted including all the four tracking times (as illustrated in Fig. 3.3). For matching based approach, the operation time is measured for all the four matching processes, which is also drawn in Fig. 3.3.

In our evaluation, 100 stereo image circular subsets captured in static scenes are used for evaluation. For a fair comparison, the number of detected feature points in each image is controlled between 400-500. Fig. 3.7 shows the results of survival rate for all the considered KLT tracker based and matching based methods. Among all the methods, KLT based STAR detector achieves the best performance (survival rate: 0.88), followed by GFTT+KLT and FAST+KLT. The BRISK detector and descriptor get the worst results. An interesting nding is that: in survival rate comparison, all the KLT based approaches perform better than feature matching based method. This fact could be explained by the repeatability in Fig. 3.5. For KLT based methods, feature points are detected once in one image and tracked to other images. While for feature matching based methods, features points are independently detected in four images. Therefore, suppose we choose SIFT detector, only 40% percent (0.8 4 ∼ 0.4) of extracted points re-appear in all the four images. This problem is not serious in KLT based methods, due to KLT tracker's high accuracy and robustness.

The results of accuracy evaluation are shown in Fig. 3.8. All the KLT tracker based methods perform better than feature matching based methods. STAR based KLT tracker achieves the best accuracy score, followed by GFTT and FAST based KLT trackers. Fig. 3.9 reveals the comparison results in terms of speed. Again, the performances of the KLT based methods are beyond feature matching based process. Matching 128-length feature descriptors of SIFT is no doubt the slowest approach. Although SURF reduce the length of descriptor to 64, ORB and BRISK step forward by substituting with a binary descriptor, they still require much more processing time than KLT tracker based method. The main reason is due to the repeated feature detection processes.

Choose the Best Approach of Circular Feature Point Association: We have compared dierent methods to achieve circular feature point association. It is obvi- ous that KLT tracker based methods are better than feature matching methods. The results are due to the characteristic of unremarkable changes between consecutive images, which makes KLT tracker more suitable than feature descriptor matching the latter has advantages when viewpoint has prominent changes. Another interesting nding is the poor performances of binary descriptors (ORB, BRISK). The improvement of matching speed for binary descriptors is maybe obtained in sacrice of accuracy, as shown in Fig. 3.8

Comparing feature detectors, we notice that STAR (CenSurE) detector performs the best in global evaluation. It is quick for computation, stable and accurate for cross-image detection and tracking. Although FAST feature detector is extremely fast for computation, it lacks robustness and accuracy for tracking, compared with STAR or GFTT. GFTT can be regarded as the second best choice, since it is just slightly worse than STAR. Therefore, in our application, we choose STAR feature detector to accurately localize feature points. KLT tracker is used for establishing point-to-point correspondences through consecutive four stereo images.

Speed Performance

Feature Detectors + Tracking/Matching method Speed (ms) 0 100 300 500 Figure 3.9: Processing speed performances for all the tested methods Technology (KIT)'s intelligent vehicle platform. The benchmark suite contains large amount of data for stereo matching, optical ow, visual odometry/SLAM, 3D object recognition. In addition, it also provides MATLAB/C++ development kits for easy access.

In our experiments, we utilize the dataset of visual odometry. The dataset consists of 11 rectied stereo video sequences with a resolution of 1241 × 376 pixels acquired with a frame rate 10f ps, in Karlsruhe. Germany. The calibration parameters and ground truth provided by a RTK-GPS are also given. The described method is tested with the dataset and several results are shown in Fig. 3.10 It can be found that, the performance of visual odometry deteriorates quickly when the cameras are turning. Also, the errors increase and accumulate along with the driving distance. More quantitative error analyzing results are shown in Fig. 3.11. Translation error analyses for these 4 sequences are given in Fig. 3.11 (a) -(h). Fig. 3.11 (i) -(j) illustrate the translation errors in average with regard to the path distance and driving speed for all the 11 sequences. It can be found from the results that, within the range of 400 meters, the translational deviation is usually around 3 -4 percentage of the total distance. However, the translation errors increase by the path distance, this is due to the errors in visual odometry coming from errors in feature points' positions, 3D triangulation and calibration parameters. In fact, long-distance and accurate visual odometry is still an open question for computer vision community. Another interesting nding is that the suitable driving speed for applying visual odometry is within 25km/h to 40km/h. The test sequences provided by KITTI benchmark are conned in a short distance. To evaluate our adopted stereo visual odometry method, we apply it using our platform SetCar, described in Sec. 1.3. The platform was driven in the city of Belfort in much longer distances (more than 1000 meters).

All the test sequences were captured by our Blumblebee XB3 stereo vision system. The stereo image pairs are well rectied after a stereo calibration process. The positions localized by the equipped RTK-GPS are regarded as ground truth and compared with the results of the described visual odometry method. Four real experimental results are shown in Fig: 3.12, 3.13, 3.14, 3.15. From the results of our experiments, it can be seen that the accuracy of the proposed visual odometry method is still unable to compare with RTK-GPS. In long distance situations, translation errors are slowly accumulated. The average translation error rate within average of 1000 meters is below 8%. The From the experiments, we can see that, compared with RTK-GPS, visual odometry still lacks of accuracy. The future works would focus on improving accuracy of visual odometry.

To improve the accuracy of feature based visual odometry, following approaches could be attempted:

• Improve the accuracy of key points locations by computing sub-pixel positions. In our experiments, the positions of feature points are in discretized coordinates. However, real world is continues. The usage of sub-pixel position of feature points maybe increase the accuracy of visual odometry.

• Model the error of key points locations. In our experiments, the weights of all feature points are the same. Modeling the error of dierent key points according to distance or other criterions maybe improve the accuracy of stereo visual odometry. In the previous chapter, we presented a method of estimating ego-motion information for a moving intelligent vehicle. However, real urban environments are complex in full of independent moving objects, such as moving pedestrians, vehicles, cyclists, etc. Therefore, the perception system of an intelligent vehicle requires an ability to analyze the moving objects when it is driven. In this chapter, we rstly present a method to detect and segment moving objects from moving stereo vision system and then, we introduce a recognition method to recognize the segmented objects based on spatial information. Vision based moving object detection is an old but still dynamic computer vision area. It provides a classication of the pixels in video sequence into either foreground (moving objects) or background. When the camera, or vision system, is static, one of the most common approach to achieve such moving object detection is background subtraction [START_REF] Wren | Pnder: Real-time tracking of the human body[END_REF][START_REF] Toyama | [END_REF]: each video frame is compared with a reference or background model, then pixels that signicantly dier from the background model are considered as parts of moving objects. The background model could be from the simplest dierence of temporally adjacent frames to complex mixture Gaussians [START_REF] Stauer | [END_REF]]. However, the problem becomes rather tricky when the vision system is moving, which is usually the case in the eld of robotics, intelligent vehicles. Similar background model is very dicult to acquire under such quickly changed environments. To cope with this problem, proposed approaches in the litterature could be roughly divided into two categories.

1. The rst category uses global motion compensation to generate a background model as utilized in motion detection in static cameras cases [Kaucic 2005]. This method suers from severe limitations in the assumption of homography transform between consecutive images. Whereas, in real moving situation, the images are changed by perspective transformation. Although several following improvements have been introduced in [Sawhney 2000[START_REF] Kang | [END_REF]]. It is still lack of accuracy and is unable to deal with complex real environment.

2. The second category is generally based on analyzing the displacement of pixels in image plane (optical ow). For example, [START_REF] Thompson | [END_REF]] incorporated dense optical ow with information of camera motion or scene structure to achieve moving object detection when the camera is also moving. [Rabe 2007a] proposed a moving object detection method based on analysis of individually tracked image points (optical ow), which provides a motion metric to indicate their motion status. In [START_REF] Dey | [END_REF]], a multiframe monocular epipolar constraint of camera motion was derived for monocular moving camera. This constraint is combined with optical ow, which captures all motion in the video across multiple frames. As an extension of optical ow from 2D image plane to 3D world, scene ow is able to detect moving objects in 3D space, as introduced in [Lenz 2011[START_REF] Wedel | [END_REF]]. This kind of methods usually involves joint estimation of ego-motion as well as objects movement. Benet of the second category is that there is no assumption for specic environment.

In this section, we present a sparse optical ow based moving object detection method. This method rstly computes sparse optical ow for detected feature points. Then, motion analysis of the sparse points is achieved by the RANSAC based visual odometry method described in Sec. 3.3: the feature points used to calculate ego-motion are assumed to be static points, while the remaining points are either from independent moving objects or noises. Based on the "suspected" feature points from moving objects, we use U-disparity image to detect real moving objects. At last, based on the detecting results in U-disparity image, the moving objects are segmented in dense disparity map. For a stereo image pair (I l , I r ), its corresponding disparity image I ∆ refers to the apparent pixel dierence or motion between I l and I r . It can be simply converted to depth map by Disparity = 1/Depth. Therefore, disparity image is usually calculated as a foundation of 3D reconstruction. Numerous algorithms have been proposed to fast and accurately compute the disparity image, for example: Block Matching [Konolige 1998], Semi-Global Block Matching (SGBM) [Hirschmueller 2008], Graph Cut [START_REF] Kolmogorov | [END_REF]], Belief Propagation [START_REF] Sun | Stereo Matching Using Belief Propagation[END_REF]]. After a balance between computational speed and quality, we choose SGBM to compute dense disparity image for our application. An example of SGBM results is shown in Fig. 4.

(b).

The gure is pseudo-colored: the color ranges from red to blue, which denotes that the disparity is from small to big. As mentioned before, our stereo image pairs are well rectied. Hence, given the disparity image, pixel-wise 3D reconstruction can be achieved by triangulation, according to Eq. 2.23. Apart from dense disparity map, U-disparity map I u (u, ∆) and V-disparity map I v (∆, v) [START_REF] Labayarde | [END_REF]] are also calculated for further scene understanding. In the eld of intelligent vehicle, U-V disparity images are practical tools for scene understanding (obstacle detection [START_REF] Wang | [END_REF]], ground plane detection [START_REF] Soquet | [END_REF]], etc.). U-V disparity images are accumulative projections of dense disparity image I ∆ to rows/columns. For example, U-disparity image is built by accumulating the pixels with the same disparity in I ∆ in each column (u-axis). Suppose a point p u = (u i , ∆ i ) in I u with intensity d u , it means that in the ith column of disparity image I ∆ , there are d u pixels with the same disparity ∆ i . Similarly, V-disparity image is obtained symmetrically. A point p v = (∆ j , v j ) in I v with intensity d v means there are d v pixels with the same disparity ∆ j in the jth row of I ∆ . Examples are shown in Fig. 4.1 (b). Actually, the U-disparity image could be viewed as a bird's view disparity image of the scene. As in the Fig. 4.1 (b), all the obstacles standing in the ground plane are mapped to a bright line in the U-disparity image. In the V-disparity image, ground plane is mapped to a quasi-line (marked as a red line). All the obstacles are projected as bright lines "grow up" from the red line. These attributes of U-V disparity images make them as convenient tools to estimate obstacles and ground plane. In our works, the U-disparity image is used for moving objects detection and segmentation, while the V-disparity image is used for estimating the ground pitch angle introduced in next chapter.

Independent Moving Object Segmentation in U-disparity Map

In most cases of driving in urban environments, the vehicle moves on a at ground plane with a certain pitch angle with respect to the vision system. Suppose in the image, there is a plane perpendicular to the ground plane. Since the ground surface is almost in the same distance to the vision system, the disparities of the pixels belonging to the ground plane are almost the same. Hence, the pixels of the ground plane in the same column are transformed to the same point in U-disparity image. This eect illustrates the reason why an object perpendicular to the ground is imaged as a bright line in U-disparity image. In Fig. 4.1, the vehicle, stone, pedestrian and pillar are all mapped as white lines in the U-disparity image. Note that the white lines in the U-disparity image are of the same column positions with the corresponding obstacles in the disparity map.

Based on the characteristic of U-disparity image, detecting moving objects is able to be achieved by detecting the white lines in U-disparity image. Recall to Chapter 3, the extracted feature points are used to compute ego-motion. Within all the associated feature points, inliers that are classied by RANSAC scheme come from the static scene, outliers correspond to independent moving objects or wrongly associated feature points. When projecting the outliers into the U-disparity image, they are always located in the bright lines. This property of U-disparity image implies that a ood-ll segmentation method [Gonzalez 1992] is feasible to detect and segment moving objects. An example is drawn in Fig. 4.2. Fig. 4.2 (a) shows the quad-associated feature points, the red points are inliers used for estimating ego-motion, the blue points are outliers coming from independent moving objects or wrongly matched feature points. Fig. 4.2 (b) shows the projection of the inliers and outliers into the U-disparity map.

Before performing ood-ll segmentation in U-disparity map, a preprocessing should be conducted to improve the performance. Restricted by the resolution of image, the moving objects too far away are usually almost the same between two consecutive frames. Moreover, the detected feature points in such area are usually not stable for analyzing. Therefore, we dene an eective range for detecting moving objects. For the used Bumblebee XB3 stereo vision system with a baseline length of 0.24m and a focal length 1007 pixels, according to Eq. 2.23, the maximum detected distance is 1007 × 0.24/1 = 241.68m. Whereas, the distance is too far away to eectively distinguish changes. In our application, we set the maximum detected distance as 30 meters. Hence, the pixels with disparities larger than 1007×0.24/30 ≈ 8 pixels are used for further analyzing.

In addition to the eective detected range, another issue should be correctly addressed. Usually, an object close to the stereo-vision system would be always captured in more pixels than an object far from the system. When transformed in U-disparity map, objects with large disparities (in small distances) are always "brighter" than objects with small disparities (in large distance), as in Fig. 4.1 (b), where the white line representing the vehicle in the left part of the U-disparity image is "darker" than the white line representing the stone block in the bottom right point, due to the distance. This phenomenon would cause diculties for further segmentation. To make the gray value distribution in U-disparity map more even, we employ an intensity adjustment as a preprocessing step. A modied sigmoid logistic function S(•) is used to adjust the intensity of U-disparity map: u are intensities of a pixel in U-disparity map before and after adjustment, respectively. r and c represent control coecients. ∆ is the row in U-disparity map. Sigmoid function is able to smoothly restrain the intensity of areas near the stereo vision system and amplify the intensity of the areas far away. Note that ∆ is bigger when an object is closer. An illustrative example is given in Fig. 4.3 (a), where three sigmoid-like functions with dierent parameters are given. With tuned parameters, an example of U-disparity map after intensity adjustment is shown in Fig. 4.3 (b), (c). We can see that the intensities of all potential objects are adjusted similar to each other, regardless of the distance.

I ′ u = I u • S(∆) = I u • r 1 + e ∆•c
Based on the corresponding feature points and adjusted U-disparity map, we segment the independent moving objects in a stereo image pair acquired at time t as follows:

1. Project all the outliers into the adjusted U-disparity map according to their disparities.

2. Take the new locations of outliers in the adjusted U-disparity map as seeds.

Then, a ood ll segmentation algorithm is performed to segment image patches with similar intensities to the seeds. After obtaining all the candidate segmentation patches, a merging process is employed to merge all the segments which are mutually overlapped.

4. Since the outliers comprise inaccurate tracked feature points appearing in static obstacles or noises, incorrect segments would lie in candidate segments.

To overcome this problem, a double-phase post-processing renement is applied. In the rst phase, each candidate segment, if it contains an inlier projection in the U-disparity map, it is rejected. In the second phase, the surviving segments are compared to stored candidate segments in previous time t -1.

If a segment has an overlapped region with a previous segment, it passes the renement and is conrmed as independent moving object in the U-disparity map. The remaining candidate segments are stored for usage in the next frame.

5. At last, conrmed segments in the U-disparity map are back projected to the dense disparity map to get independent moving objects in the dense disparity map. To better evaluate the proposed method, we quantitatively analyze the moving detecting results by comparing them with ground truth. In this experiment, ve test sequences (each of them is composed of more than 100 frames) are used for evaluation. We compute the Type I errors (true positive rate and false positive rate) for each sequence, the results are as follows: From Tab. 4.1, we could see that the proposed framework performs well when there are moving vehicles in the scenario. However, the performance degenerates when pedestrians appear. The major reason is that a moving pedestrian is slower than a moving vehicle. Hence, several detected feature points within a moving pedestrian are classied as static. This problem would cause undetected moving pedestrians. In the previous sections, ego-motion of the vehicle itself and independent moving objects are investigated. In this section, we propose a recognition method to classify moving objects (such as pedestrian, vehicle, cyclist) by sparse spatial information.

Introduction

Vision-based vehicular perception systems usually consist of two distinct levels of functions which rely on two distinct levels of information, respectively. For the rst level, spatial information, either dense or sparse, is always utilized to deal with the problems of motion analysis, including ego-motion estimation, moving object detection and tracking [Nedevschi 2007[START_REF] Rodriguez | [END_REF]] [Rabe 2007b]. For the second level, some complex appearance features are usually exploited in object recognition and semantic segmentation [Dalal 2005[START_REF] Viola | [END_REF][START_REF] Leibe | [END_REF]].

In the last decades, one trend of vision-based perception systems for intelligent vehicles is breaking through the boundary of the two levels. In [Leibe 2008a], the authors integrated both structure information and complex appearance features (implicit shape model (ISM) detector [Leibe 2008b]) for object detection and tracking.

In [Hoiem 2008], the authors improved vehicle, pedestrian detection by considering 3D perspective information. In [Brostow 2008], semantic segmentation of trac scenarios is completed by only structure and motion information. Classication of particular types of objects is a dynamic research eld. A common scheme of solution contains two steps: At rst, extracting and processing certain appearance features (as SIFT [Lowe 1999], HOG [Dalal 2005], Haarlike features [START_REF] Viola | [END_REF]]) from labeled dataset. Second, a supervised machine learning method (as SVM [Burges 1998], AdaBoost [START_REF] Freund | [END_REF]], Random Forest [Breiman 2001], Gradient Boosting [Friedman 2000]) is then applied to train a good classier for further prediction. However, since the previous approaches heavily rely on appearance features, several researchers tried to introduce geometric information into this eld. In [Brostow 2008], the author improves semantic segmentation and recognition by taking advantages of 3D geometric features. In [Hoiem 2008], vehicle and pedestrian detection results are improved with the help of perspective information. In [START_REF] Zhu | [END_REF]], a range image generated by a laser range nder is segmented and classied by only depth information.

Our proposed method is inspired from [Brostow 2008] and [START_REF] Zhu | [END_REF]]. Our objective is attempting to classify moving objects in urban environments solely by spatial information. The spatial information is generated from segmented pixels in stereo images. Afterward, several basic spatial features as well as advanced spatial features (based on Kernel Principal Component Analysis (KPCA) [Schölkopf 1998]) are extracted. Finally, a classier which has been trained o-line, is used to predict the categories of the moving objects.

Spatial Feature Extraction and Classication

After the previous steps, a segment of moving object is extracted as a set of pixels. By triangulation, it can be represented as a cloud of 3D points: C = {P 1 , ..., P M }. In this section, several spatial features are extracted and trained to recognize common moving objects in urban environments.

Spatial Feature Extraction

The authors in [Brostow 2008] proposed 5 simple structure features of every 3D point for semantic segmentation. For a cluster of 3D points, [START_REF] Zhu | [END_REF]] utilized 5 simple structure features and 2 statistical features. Although improvements are achieved by these methods, they still leave large space to work in spatial features. In fact, the limitations of their spatial features arise from the intuitive representations of points in Euclidean space. The statistical normal features in [START_REF] Zhu | [END_REF]] is in fact the main directions of a 3D points cluster which is equivalent to doing Principal Component Analysis (PCA) for the points cluster. Traditional PCA is a powerful method for extracting structures from high-dimensional data. However, data of a 3D points cluster is low-dimensional (Euclidean space).

(a) Two linearly inseparable co-centric point clusters (b) The projections of the two clusters in a transformed higher dimensional space using the rst two principal components With the aid of kernel methods, more and advanced spatial features can be extracted by transforming the feature space from 3D Euclidean space into a higher dimensional space. The reason stems from that two linearly inseparable clusters can sometimes become linearly separable in a transformed higher space, as shown in Fig. 4.6. The left image shows two clusters which can not be divided by a linear method, while after transforming them into a higher space (by Gaussian kernel), the two clusters are easily to be separated from their projections into the rst two principal components. In fact, if N points cannot be linearly separated in d < N dimensional space, they always can be almost linearly separated in d ′ ≥ N dimensional space. For instance, given N points P i , i = 1, ..., N in d dimension, if we map them into an N -dimensional space with by a mapping Φ(•):

Φ(P i ) = δ ij (4.2)
where Φ : R d → R N and δ ij is the Kronecker delta. It is not dicult to construct a hyperplane that divides the points into arbitrary clusters (see in Fig. 4.6 (b)).

Inspired by this phenomenon, we utilize a Kernel PCA (KPCA) based method [Bo 2011] to improve the features in [START_REF] Zhu | [END_REF]]. The method rstly transforms the 3D points cluster into a higher dimensional space by KPCA. Then, the histogram of top L eigenvalues of kernel matrix is calculated and taken as advanced spatial features. KPCA is an extension of PCA using kernel methods to map data from low dimension into high dimension. For a points cluster C = {P 1 , ..., P M }, where

P i = (X i , Y i , Z i ) T .
The steps of KPCA are as follows:

1. Building the M × M kernel matrix K, where:

K(i, j) = (Φ(P i ) • Φ(P j )) = Φ(P i ) T Φ(P j ) (4.3)
In practice, the commonly used kernel functions are: polynomial kernel: k(x, y) = (x • y + a) c Gaussian kernel: k(x, y) = exp(-γ∥x -y∥ 2 ), exponential kernel: k(x, y) = exp(-γ∥x -y∥).

2. Centering the kernel matrix by:

K = K -1 M K -K1 M + 1 M K1 M (4.4)
where

1 M is a M × M matrix with (1 M ) ij = 1/M .
3. SVD decomposition of the centered kernel matrix:

K = USV * (4.5)
4. In the matrix S, the singular values of K are arranged in the diagonal with decreasing order. The top L singular values together with their corresponding column vectors in S represent the main directions of a transformed points cluster. Hence, two kinds of objects in dierent shapes would have dierent distribution of these main directions, as shown in Fig. 4.7. In practice, the top L singular values λ i , i = 1, ..., L are selected and normalized as:

f i = λ i / L ∑ i λ i (4.6)
In Table . 4.2, all the features extracted from the cluster in practice are listed.

Classication

The objective of classication aims to recognize the common moving objects in urban environments, including vehicles, pedestrians and cyclists .... Based on the o-theshelf machine learning methods, together with the spatial features extracted in the previous step, we adopt two popular machine learning methods (Random Forest [Breiman 2001] and Gradient Boosting Trees [Friedman 2000] [Friedman 1999]) to train classiers based on the extracted spatial features.

Random Forest: Random forest is an ensemble learning framework for classication (or regression) that operates by constructing a multitude of decision trees at training time. The rst formal introduction of random forest was given in [Breiman 2001], which describes a method of building a forest of uncorrelated trees using classication and regression trees(CART) [START_REF] Breiman | [END_REF]] like decision trees, combined with randomized node optimization and bagging. In his paper, each tree is constructed using the following steps:

• Let T be the number of training cases, and the number of variables in the classier be V .

• v denotes as the number of input variables to be used to determine the decision at a node of the tree; v ≪ V .

• Choose a training set for this tree by choosing n times with replacement from all T available training cases (e.g. take a bootstrap sample). Use the rest of the cases to estimate the error of the tree, by predicting their classes.

• For each node of the tree, randomly choose m variables on which to base the Recognition decision at that node. Calculate the best split based on these m variables in the training set.

• Each tree is fully grown and not pruned (as may be done in constructing a normal tree classier).

For prediction, a new sample is pushed down the tree. It is assigned the label of the training sample in the terminal node it ends up in. This procedure is iterated over all trees in the ensemble, and the mode vote of all trees is reported as the random forest prediction. Compute pseudo-residuals:

r im = -[ ∂L(y i ,F (x i )) ∂F (x i )
] F (x)=F m-1 (x) f or i = 1, ..., n.

5:

Fit a base learner h m (x) to pseudo-residuals.

6:

Compute multiplier γ m by solving the following one dimensional optimization problem:

γ m = arg min γ n ∑ i=1 L(y i , F m-1 (x i ) + γh m (x i )) (4.8) . 7:
Update the model:

F m (x) = F m-1 (x) + γ m h m (x) 8: end for 9: Output F M (x).
Gradient boosting trees is a special case of gradient boosting by using x-sized decision trees (especially CART trees) as base learners.

4.2.3

Experimental Results

Dataset

We tested the proposed recognition method both on KITTI dataset and our dataset.

All the stereo images are rectied and the intrinsic and extrinsic parameters of the cameras are also provided. The segmentation results from Sec. 4.1 are used for further recognition. For each segment, we reconstruct its pixels' 3D positions by triangulation. The spatial features are extracted according to Tab. 4.2 and sent to trained classiers for recognition. Three commonly used kernel functions (Polynomial, Gaussian and Exponential Kernel) were adopted. where c = 2, a = 1, γ = 0.5 (see Sec. 4.2.2). We trained two types of classiers (Random Forest, Gradient Boosting Trees) using the training dataset. The parameters of the two kinds of classiers (as the depth of a decision tree, the maximum number of trees, the loss function type ...) are optimized by cross validation between the training and test data. Moreover, we also implemented the proposed features proposed in [START_REF] Zhu | [END_REF]] for a comparision. Since the detailed implementation of [START_REF] Zhu | [END_REF]] is not given, we simulate his method as well as we can. All the results are shown in Table . 4.4, which illustrates:

• In theory, Random forest is an example of bagging approach, less prone to overt. Gradient boosted trees represent the other one. Both are very successful in many applications. While in our case, the Gradient Tree Boosting classier performs better than Random Forest classier in all the cases in our dataset.

• The classication results based on the three kinds of kernels are almost the same. However, the Gaussian kernel performs slightly better than the other two kernels.

• Compared with Zhu's method [START_REF] Zhu | [END_REF]], the proposed kernel shape descriptors improve the classication results thanks to transforming the features into a higher dimensional space. In average, the precision based on the proposed kernel shape descriptors increases by around 10% than in [START_REF] Zhu | [END_REF]].

• The Confusion matrices reveal that pedestrians and cyclists samples are apt to be confused. This is because they are similar in both size and shape.

The experimental results for the moving object recognition are shown in Fig. In this chapter, we proposed stereo vision based methods to detect, segment and recognize independent moving objects, for a moving intelligent vehicle. These methods are based on the ego-motion estimation results in Chap. 3.

In the rst part, we presented a stereo vision based method to detect and segment independent moving objects from a moving intelligent vehicle. The method relies on the "outliers" generated in the step of visual odometry. The "outliers" mainly consist of matching noises and feature points from moving objects. Next, they are projected to U-disparity image as seeds for moving object detection and segmentation. This method is tested on our datasets and the experimental results are analyzed. However, the performance of this method degenerates when trying to detect moving pedestrians. The main reason of this problem is that moving pedestrians usually move too slow to eectively detect their motion between two consecutive frames. Future works would be around this problem. One possible solution to this problem is using pedestrian detection methods detect pedestrians at rst. Then, if there are any detected moving feature points in detected pedestrian region, it could be classied as a moving pedestrian.

In the second part, a recognition approach based on learning spatial information is proposed. We tried to classify the segmented moving objects acquired by the previous step. In this method, a kernel PCA algorithm is used to select the spatial features in high dimensional space. Next, we trained two kinds of classiers: Random Forest and Gradient Boosting Trees to recognize three kinds of moving objects: Pedestrian, Vehicle and Motor/Bike. The experimental results showed good recognition rates, especially for the Gaussian kernel based KPCA features combining with Gradient Boosting Trees. Also, we show that with the boosted kernel PCA features, our method outperforms an another method by around 10% in terms of precision. In future works, we consider combining spatial features with appearance features to reach a better recognition results. For example, the histogram of oriented gradients (HoG) or local binary patterns (LBP) could be introduced to improve the recognition performance. Introduction As we introduced in the rst chapter, we plan to use stereoscopic system and 2D lidar to analyze the environment for an intelligent vehicle. However, despite many impressive progresses in both computer vision, and range sensors in past decades, no technique based on single kind of sensors can be reliably and robustly performed in a Lidar real complex environments. In fact, visual sensors are always limited in narrow eld of views (FOV) and require huge computational cost. On the other hand, providing only depth information, the range sensors are not enough for more complex tasks in deep level. Hence, achieving eective and robust environments understanding requires the fusion of dierent kinds of sensors. Indeed, fusing visual and range data has attracted more and more attentions for intelligent vehicles. In [START_REF] Mahlisch | [END_REF][START_REF] Fox | A spatio-temporal probabilistic model for multisensor object recognition[END_REF][START_REF] Hwang | [END_REF], various information fusion methods are proposed.

The problem of fusing data between visual cameras and LIDARs arises from the heterogeneity of the sensors. Lidar data is depth information in a certain scan plan, while the data of a camera can be viewed as a projection of a 3D world into an image plane. Thus, how to bridge the gap between a visual camera and a LIDAR has to be solved before performing data fusion. An initial solution is to estimate a rigid geometrical transformation (a 3 × 3 rotation matrix and a 3 × 1 translation vector) between lidar coordinate system and stereo vision coordinate system which is usually called extrinsic calibration between lidar and stereo vision system. In general, it is implemented as the rst step for further information integration of the two sensors, as illustrated in [START_REF] Perrollaz | [END_REF][START_REF] Cristiano | [END_REF], Huang 2009b].

Related works

In literature, most proposed methods for extrinsic calibration between a lidar and a camera can be roughly divided into three categories:

1. Methods based on auxiliary sensors. In [START_REF] Aliakbarpour | An Ecient Algorithm for Extrinsic Calibration between a 3D Laser Range Finder and a Stereo Camera for Surveillance[END_REF]], with the aid of an IMU (Inertial Measurement Unit) and a visible laser point, the Euclidean 3-dimensional transformation between a lidar and a camera is acquired by rotating the scan plane, where the angular information is estimated by the inertial sensor. Nunez et al. [START_REF] Nunez | [END_REF]] also use an IMU to construct geometric constraints between a lidar and a camera system for calculating the transformation. However, being equipped with an IMU is a prerequisite for this kind of approaches.

2. Methods based on specially designed calibration boards. Li et al. [Li 2007] design a right-angled triangular checkerboard and employ the invisible intersection points of the lidar's scan plane with the edges of the checkerboard to set up constraint equations. Rodriguez et al. [START_REF] Rodriguez | [END_REF] adopt a calibration board with two concentric circles where the inner circle is hollowed. The authors use the correspondence of the center of the inner circle between several dierent poses to achieve extrinsic calibration. By making use of special features in the designed calibration board, such approaches enable one to acquire the geometric relationship between a lidar and a camera. However, specially designed calibration boards impede themself for more widespread uses.

Methods based on common calibration chessboards.

Zhang and Pless [Zhang 2004] rstly proposed a method utilizing an ordinary calibration chessboard which is widely used in camera calibration. This method is based on the fact that lidar measurements on the calibration chessboard satisfy certain geometric constraints with the plane of the chessboard. Following Zhang's idea [Zhang 2004], Huang and Barth [Huang 2009a] utilize similar geometric constraints and generalize them into a situation of multi-layer lidar. Vasconcelos et al [Vasconcelos 2012] improve Zhang's method [Zhang 2004] by reducing the minimum number of various poses from 6 to 5. This kind of methods can handle with intrinsic calibration for cameras and extrinsic calibration between a lidar and a camera at the same time. Also, these methods require no extra equipments or specially designed calibration boards. Thus, in our point of view, they are more practical in real applications.

According to the above analysis, we prefer the third approach for its simplicity. Since stereovision can provide a complete three-dimensional view of a scene, its application for intelligent vehicles is more and more comprehensive nowadays. However, all the kinds of methods mentioned above are applied between a monocular camera and a 2D lidar. Calibrating a stereo vision system with a lidar can be performed by calibrating each camera with the lidar separately, as all existing methods perfom, or more reasonably, by calibrating the considered stereovision system with the lidar. Based on 3D reconstruction of a common calibration board, we propose a novel calibration approach for obtaining the relative position and orientation between a stereo vision system and a lidar. It is worthy to mention that the advantages of our approach are threefold: Compared with the rst and second categories of methods, our method is more convenient since it uses one standard chessboard; Compared with the third category, we will show that owing to the consideration of the whole stereovision system, the proposed extrinsic calibration algorithm performs more precisely. Moreover, the introduction of Mahalanobis distance contraints and the consideration of sensor noise models make the calibration results more accurate and robust.

Problem Formulation

A geometric model of the multisensor system together with the 2D chessboard is given in Fig.

= R 3 lef t . R 3
right is linked to the stereoscopic coordinate system R 3 stereo by a rigid transformation composed of a rotation matrix R l r and a translation vector T l r . Suppose a point P in the right camera coordinates as: P r = (X r Y r Z r ) T , then, its corresponding coordinates in the stereoscopic system R 3 stereo are: In a similar way, a point in the stereoscopic coordinate system can be mapped into the right camera coordinate system as Eq. ( 5.1), by substituting R l r and T l r as:

P stereo = P l = R l r • P r + T l r ( 5
R r l = R l r T T r l = -R l r T • T l r (5.2)
Let (Φ l , ∆ l ) and (Φ r , ∆ r ) be the 3D rigid transformation from the LIDAR coordinate system R 3 lidar to the left and right camera coordinate systems R 3 lef t and R 3 right respectively, where Φ l and Φ r are orthogonal rotation matrixes, ∆ l and ∆ r are translation vectors. Suppose a 3D point P observed by both the LIDAR and the stereoscopic system, and denoted as P lidar = (X lidar Y lidar Z lidar ) T in R 3 lidar , P l and P r in the left and right camera coordinate systems. The three coordinates are connected by:

P l = Φ l • P lidar + ∆ l P r = Φ r • P lidar + ∆ r P l = R l r • P r + T l r (5.3) 5.1.3 Performance Evaluation
In practice, it is dicult to obtain the ground truth of the real extrinsic parameters between a camera and a lidar. In [Zhang 2004], [Huang 2009b], the extrinsic calibration methods are evaluated just by intuition (whether the projected laser points in images are reasonable or not). In this paper, we dene an indicator ε as a precision measure of the multi-sensor calibration system:

ε 2 = 1 L L ∑ i ∥R l r • (Φ r P i lidar + ∆ r ) + T l r -(Φ l P i lidar + ∆ l )∥ 2
(5.4) where L refers to the number of lidar points on the chessboard used in the extrinsic calibration. ε is not a direct error with respect to the ground truth. Indeed, it can be viewed as a precision measure evaluating the conformity of the results with respect to the geometric relationships within the multi-sensor system. We hope ε to be as small as possible.

3D Plane Reconstruction based Extrinsic Calibration

A owchart in Fig. 5.2 illustrates all the steps of our method. It requires putting a calibration chessboard in front of all the sensors at dierent positions and orientations and make sure that all the sensors can detect it. For each pose, acquired a Lidar images are undistorted and rectied as in Sec. 2.2.3.2. Afterwards, we automatically extract the corner points of the chessboard in the images of the two cameras, and triangulate their 3D positions in the stereoscopic coordinate system. Then, from the reconstructed points, a PCA based method is used to estimate the 3D plane, which best ts the reconstructed points of the chessboard in the stereoscopic coordinate system. Finally, considering the geometrical constraints, a non-linear 3D minimization is carried out in order to calculate the extrinsic parameters. Our stereoscopic system delivers rectied stereo image pairs automatically, all the intrinsic/extrinsic parameters of the stereo-rig are calculated by [Zhang 2000]'s method (implemented in [START_REF] Bouguet | Camera Calibration Toolbox for Matlab[END_REF]). Therefore, the intrinsic/extrinsic parameters K l , K r , R l r , T l r , are known and K l = K r .

Corner Points Triangulation

In Sec. 2.2.4, a two-view triangulation method in ideal situation is given. However, in general cases, the two cameras can be rotated with respect to each other. Furthermore, due to the inevitable imaging noises, the two rays O l P and O r P are not guaranteed to be intersected. The process of 3D triangulation becomes complicate in such cases. A common method (middle point method) [Beardsley 1994] is to use the middle point of the perpendicular to the two rays, as drawn in Fig. 5.3 (a). Beyond a simple middle point method, [START_REF] Hartley | [END_REF]] proposes an optimal triangulation method to minimize projection errors of corresponding points under epipolar constraint. In [Kanatani 2008], a method is proposed to correct the positions of corresponding to achieve optimized triangulation results. Fig. 5.3 (b) illustrates the case of optimal triangulation. More details about 3D triangulation in general cases could be found in [Hartley 2004].

In 3D computer vision practices, rectifying stereo image pairs is a necessary preprocessing step to simplify the succeeding 3D reconstruction. In our approach, input stereo images are all undistored and rectied. Afterward, we apply a sub-pixel corner point extraction [START_REF] Bouguet | Camera Calibration Toolbox for Matlab[END_REF]] to get sets of corresponding pairs of corner points p l = (u l , v l ) and p r = (u r , v r ). 3D reconstruction is performed according to Eq. 2.23. Let P l be the 3D point coordinates in the left camera coordinate system, P r be its corresponding coordinates in the right camera coordinate system.

By triangulating all the corner point pairs extracted from the two views, the corresponding 3D coordinates attached to the left and right camera coordinate systems are acquired. Although errors exist in the results of triangulation, the 3D points cloud is nearly at. The next step is to estimate a 3D plane that best ts all the reconstructed 3D points.

3D Plane Estimation

A plane in 3D Euclidean space can be dened by a known point on the plane and a normal vector perpendicular to its surface. In homogeneous coordinate system, the where

⃗ n = [a, b, c] is a 1 × 3 normal vector ( √ a 2 + b 2 + c 2 = 1), d = -⃗ n • [x 0 y 0 z 0 ] T , [x 0 y 0 z 0 ]
is a known point on the plane. d represents the perpendicular distance from the origin to the plane. More simple, Eq. 5.5 can be written as:

ax + by + cz + d = 0 (5.6)
In mathematics, given a set of 3D points S = {P 1 , ..., P i , ..., P N |i = 1, ..., N }, evaluating its best tting plane is a classic linear regression problem. As introduced in Sec. 2.3.1, linear least squares methods can be applied to solve this problem. Ordinary Least Squares Plane Fitting: As described in Sec. 2.3.1.1, ordinary least squares (OLS) assumes that only the predicted variables have errors. In 3D plane tting, the problem becomes to nding plane coecients to minimize a summation of residuals:

N ∑ i=1 |z i -(a ′ x i + b ′ y i + d ′ )| 2 (5.7) where a ′ = -a/c, b ′ = -b/c, d ′ = -d/c.
For convenience, the plane equation represents in Eq. 5.6 is changed as

z = a ′ x + b ′ y + d ′ (5.8) a Lidar
According to Eq. 2.32, let P = [x, y, 1] T be a N × 3 observational matrix, z be a N × 1 predicted vector. The OLS solution of 3D plane tting is:

  a ′ b ′ d ′   = (P T P) -1 P T z
(5.9) Supposing the predicted variable z is contaminated by independent and identically distributed Gaussian noises with variance σ, the covariance of the plane coecients is given by Eq. 2.34 as:

Q = (P T P) -1 σ 2
(5.10) Total Least Squares Plane Fitting: Recall to Sec. 2.3.1.3, total least squares extends the assumption of ordinary least squares by considering errors in all variables. According to Eq. 2.49, TLS computes a right singular vector corresponding to the minimum singular value of an augmented matrix. The TLS model for 3D plane tting is:

z + n z = (a ′ , b ′ , d ′ ) • (   x y 1   +   n x n y 0   ) (5.11)
where n x , n y , n z are the errors in all the observational data (X, Y, Z). Recall to Eq. 2.48, the TLS 3D plane tting problem is:

minimize nx,ny,nz ∥[n x , n y |n z ]∥ F subject to z + n z = (a ′ , b ′ , d ′ ) • (   x y 1   +   n x n y 0   )
(5.12) Let C = [x, y, 1, z] be a N × 4 augmented matrix comprising all the observational data. Then, we perform a singular value decomposition to C:

C = UΣV T (5.13)
where U and V are 4 × 4 real unitary matrices. Σ is a 4 × 4 diagonal matrix with eigenvalues on the diagonal in decreasing order:

λ 1 ≥ λ 2 ≥ λ 3 ≥ 0. Let v i be the ith column of V.
Then, v 4 is the eigenvector corresponding to the minimum eigenvalue. According to Eq. 2.49, the TLS solution for 3D plane tting is:

  a ′ b ′ d ′   = - t α (5.14)
where t is the rst three elements of v 4 , α is the fourth element of v 4 . Assuming the errors in all observational data (X, Y, Z) are independent identical distributed zero-mean Gaussian noises: n x , n y , n z ∼ N (0, σ), the covariance matrix of TLS estimation results is approximated according to Eq .2.50:

Q = 1 γ σ 2 ∥(a ′ , b ′ , d ′ )∥ 2 (M T M) -1
( 5.15) where M = [x, y, 1] is a N × 3 matrix. In fact, the TLS method in 3D plane tting equals to principal Component Analysis (PCA) [START_REF] Wall | Singular Value Decomposition and Principal Component Analysis[END_REF]].

Applying the method described above, 3D plane estimation is performed for various positions and orientations of the calibration board. Fig. 5.4 illustrates an example of six 3D planes estimated according to six dierent sets of reconstructed 3D points. The black frames represent the estimated planes and the colored points within the frames are the reconstructed 3D corner points. a Lidar

Automatic Lidar Measurements Extraction

In [Zhang 2004] and [Huang 2009b], the ways of extracting lidar points on the calibration board are not mentioned. Here, we apply an automatic extraction approach by dierencing the acquired measurements and background data in a static environment. The background lidar measurements are acquired before starting calibration.

Then, the measurements acquired during the calibration process are compared with the background by a simple subtraction operation to get the lidar points on the chessboard. According to [START_REF] Ye | Characterization of a 2-D Laser Scanner for Mobile Robot Obstacle Negotiation[END_REF]], when a laser ray hits on the edge of the chessboard, the measurement would be the average distance between the chessboard and the background. Consequently, the lidar points in the edge of the chessboard which have unreasonable measurements should be discarded. The nal results are shown in Fig. 5.5. The steps presented above allow us to estimate 3D planes of the calibration board with various poses (j = 1, ..., M ). In ideal situations, every lidar measurement on the calibration board should satisfy the plane equation of the calibration board:

⃗ N • Pstereo = ⃗ N • H • Plidar = 0
(5. 16)

where H = [ Φ l ∆ l 0 1 ]
denotes 3D rigid transformation which we want to estimate. , c, d] is the estimated 3D plane of the chessboard. However, in real experiments, Eq. 5.16 is rarely satised. Given an estimated rigid transformation H, all the residuals can be represented as:

⃗ N = [⃗ n, d] = [a, b
e = { ⃗ N 1 H P 1,1 lidar . . . ⃗ N 1 H P L 1 ,1 lidar . . . ⃗ N M H P 1,M lidar . . . ⃗ N M H P L M ,M lidar }
(5.17)

where P k,j lidar |(k = 1, ..., L j , j = 1, ...M ) is the homogeneous coordinate of the kth lidar point on the chessboard at the jth pose. ⃗ N j |(j = 1, ...M ) is the 3D plane of the chessboard put at the jth pose.

Although directly minimizing a summation of squared errors in Eq.( 5.17) makes sense, it doesn't take into account enough the errors produced either from the 3D plane equations estimated by the stereovision system or from the lidar itself. Considering the error of 3D triangulation grows quadratically, this could make results worse when the calibration board is far away from the stereovision system. To make the method more robust against errors, we propose to minimize a summation of squared Mahalanobis distance, rather than Euclidean distance. Compared with Euclidean metric, the Mahalanobis distance benets can be viewed as follows: rstly, it corrects the scale factor. In our work, the errors come from both stereoscopic system and lidar. Euclidean distance is sensitive to the scale of various kinds of variables involved. For Mahalanobis distance, thanks to the measure of covariance, the problem of scale factor is resolved. Secondly, Mahalanobis distance corrects the correlation between variables. This advantage is also a consequence of calculating covariance. Assuming e ∽ N (0, C), then the optimization of extrinsic parameters based on Mahalanobis distance becomes:

H(Φ l , ∆ l ) = arg min Φ l ,∆ l {e T C -1 e} (5.18)
where C is the covariance matrix of e. Assuming the elements of e in Eq.( 5.17) are a Lidar mutually independent, then C is a diagonal matrix as:

C =              C 1,1 . . . 0 C L 1 ,1 . . . C 1,M 0 . . . C L M ,M             
( 5.19) where C k,j = var( ⃗ N j H P k,j lidar ). Meanwhile, assuming the plane estimation ⃗ N and measurement P lidar are two independent random vectors with multinormal distribution N ( ⃗ N, Σ ⃗ N ) and N ( Plidar , Σ P lidar ). A closed-form solution for every element of C is complex. Noticing that we don't know H before computation. Meanwhile, H plays a role as scale factor when expanding C k,j . Here, we approximate C by ignoring the middle term H:

C k,j ≈ var( ⃗ N j P k,j lidar ) = ⃗ N T j Σ P k,j lidar ⃗ N j + P k,j T lidar Σ ⃗ N j P k,j lidar + T r(Σ P k,j lidar Σ ⃗ N j ) (5.20)
where T r( * ) is the trace of a matrix. Hence, the covariance matrix can be approximated with the distributions of N ( ⃗ N, Σ ⃗ N ) and N ( Plidar , Σ P lidar ). The covariance matrix of lidar measurement is inferred in Sec. 2.1.2 and given in Eq. 2.6. For the covariance matrix of estimated 3D plane coecients, Eq. 5.15 gives a covariance matrix for the coecients (a ′ , b ′ , d ′ ). Since in our work, the plane coecient c is always set to -1, corresponding covariance matrix of plane coecients is: (5.21) where Q ij refers to the (i, j)th element in Eq. 5.15.

Σ ⃗ N =     Q 11 Q 12 0 Q 13 Q 21 Q 22 0 Q 23 0 0 0 0 Q 31 Q 32 0 Q 33    
After obtaining the covariance matrix C in Eq. 5.19, we adopt the Levenberg-Marquardt algorithm to nd the optimized H, which minimizing the summation of residuals as in Eq. 5.18.

Summary of the Calibration Procedure

The proposed calibration procedure is summarized below in the form of Algorithm 5.

Algorithm 5 3D Plane Reconstruction Based Extrinsic Calibration 1: for j = 1 to M (number of poses) do 2:

Reconstruct their 3D position in the stereoscopic coordinate system by optimal triangulation, as presented in section 5.2.1

3:

Estimate the best tting plane [ ⃗ N j , d j ], as stated in section 5.2.2 4:

Automatically select the points {P k,j lidar |k = 1, 2, ..., L j } on the calibration board from lidar measurements 5:

For every selected lidar point and the estimated plane, compute their covariance matrixes according to Eq. (2.6) and Eq. ( 5.21), respectively 6: end for 7: Given a rst guess, [Φ 0 , ∆ 0 ], use LM-algorithm to reach a convergence with [Φ, ∆] minimizing the objective function of Eq. ( 5.18) Gaussian noise with zero mean and standard deviation σ (from 0.5 to 6 pixels) is added to the virtual corner points. Based on those "noised" corner points, the proposed method is compared with Zhang's method [Zhang 2004]. The errors are computed as: The lidar points are "polluted" by zero-mean Gaussian noises of range and angle measurements n r and n θ with variances δ r and δ θ , respectively. δ θ is set to 0.0001 degree while δ r varies from 0 to 50 (mm). The simulation results are shown in Fig. 5.8. One can see that our method makes the calibration results more robust, when compared with Zhang's method [Zhang 2004]. We randomly generate 22 dierent 3D calibration board planes 10 times. By gradually increasing the number of various poses of the chessboard from 2 to 22 in each run, we record the results and calculate the errors. We also test Zhang's method [Zhang 2004] for comparison. We run the tests 10 times with various sets of chessboard poses. The results are averaged, and shown in Fig. 5.9. The error bars represent the variance of 10 dierent tests. Good results can be obtained when the number of poses is larger than 6.

ε Φ l = ∥Φ l -Φ ′ l ∥, ε ∆ l = ∥∆ l -∆ ′ l ∥ (

t the rst guess in LM algorithm

Although Levenberg-Marquardt (LM) algorithm is very robust, it is required to evaluate the inuence of the rst guess on the nal results. In this experiment, since it is impossible to exhaustively try all the random rst guesses, we test all the factors in a rst guess independently. For the initial rotation guesses, we start from a 3 × 3 identity matrix, then independently change the yaw, pitch and roll angles [wik 2012] from 0 • to 90 • with an interval of 5 • as:

R guess =   1 0 0 0 1 0 0 0 1   • R x (γ) • R y (α) • R z (β) (5.23)
where γ, α, β are the yaw, pitch and roll angles. For the initial translation guesses, we start from [0, 0, 0] T (in meter), then gradually increase the X, Y, Z values independently to [10, 10, 10] T (in meter) with an interval of 0.5 meter. Results of errors according to various rst guess values are shown in Fig. 5.10. It is shown that LM algorithm is quite robust. Indeed, even for a rst guess that is far from the expected value, it always produces solutions close to the ground truth. The calibration process is performed in an indoor environment at rst, then, the results are veried in outdoor experiments. The calibration pattern we made is of 19×12 squares with equal size 50mm×50mm. The chessboard is detected by the two sensors meanwhile with 20 various poses. We both apply our method and Zhang's approach [Zhang 2004] for comparison. The projection of the lidar measurements into the images of the stereoscopic system is shown in Fig. 5.11. It seems that .11: Comparison of the proposed method with Zhang's method [Zhang 2004]. In (c) and (d), the blue points (by the proposed method) are almost at the same height, while the red points (by Zhang's method [Zhang 2004]) are not.

both of the results are quite reasonable. While from the partial enlargement views, our method gets more stable and precise results. Indeed, blue points (calculated by our method) in the pair of images have almost the same height, while the red ones (calculated by Zhang's method [Zhang 2004]) have obvious deviation between the two images. The estimated extrinsic parameters using our method and Zhang's method are shown in Tab. 5.1. It shows that our algorithm performs more precisely than Zhang's method [Zhang 2004]. The measure ε in Eq. 5.4 explains the deviation of red points in Fig. 5.11. The proposed method Zhang's method [Zhang 2004] Φ l [0.1474, -0.0251, 0.0083] [0.1520, -0.0284, 0.0045] ∆ l [0.1195, -1.4543, -1.2328] [0.1245, -1.4566, -1.2297 More experiments are used to test the calibration method in outdoor environments, as shown in Fig. 5.12 and Fig. 5.13. The experiments are performed in the city of Belfort, France. Left columns of Fig. 5.12 and 5.13 show six typical scenes in urban area captured by the left camera, as well as the corresponding lidar measurements within its view. Right columns of Fig. 5.12 and 5.13 show lidar measurements (the 2D lidar mounted on the bottom of our platform). The red points are within the view of the left camera, and then are projected into corresponding images as shown in the left column.

We also apply the same algorithm for the 2D lidar mounted on the top of our platform. Fig. 5.14 (b),(c) show lidar measurements from the bottom and top. Fig. 5.14 (a) shows the projections of lidar points into the left camera's eld of view by the output of the proposed algorithm. Another similar example is shown in Fig . 5.15. The lidar points within the camera's eld of view are marked by red color. We can see that the projections are quite reasonable, since the geometric characteristics of the scenes detected by the camera coincide closely with lidar data.

Conclusion And Future Work

In this chapter, we presented a novel extrinsic calibration method for integrating a lidar and a binocular stereoscopic system. Since the geometric relationship between the two cameras in the stereoscopic system is considered, real data experiments showed that our method is more exact and precise than Zhang's method which is widely used. Meanwhile, the introduction of sensor noise models and non-linear Mahalanobis distance optimization makes our method much more robust than Zhang's method. This is veried by computer simulation results. The described method can be applied for any multi-sensor fusion system consisting of multiple cameras and multiple lidars. Future work will be mainly focused on the calibration of a stereovision system with multi-layer lidars. Also, we plan to develop an improved self-calibration method.

In the eld of intelligent vehicles, many tasks, such as localization, collision avoidance, path planning, are usually performed based on well represented maps [Leonard 2008[START_REF] Nguyen | [END_REF]. Mapping is the process of transforming sensor measurement to an image of the environment. In the early days of robotic research, many tasks, such as localization/navigation, depend on the map being dened as a nite sized set of landmarks. However, physical sensors do not usually detect landmarks unambiguously. In order to use landmark based algorithms, sensor readings must be pre-processed in a separate step to convert the raw sensor data into a set of detected landmarks, such as in [Leonard 1991]. The additional step introduces more error into any algorithm, as well as discarding much of the sensor information which does not contain any landmark.

One of the primary drawbacks of landmark based maps is data association problem. Because raw sensor data are not labelled with correct landmarks detected, the Figure 6.1: Example of occupancy grid map [Miroslav 1999] sensor processing must somehow determine exactly which landmarks were observed. If mistakes are made, the localization and mapping algorithms, which depend on sensor data, will fail. In order to compensate the data association problem, many localization and SLAM (Simultaneous Localization And Mapping) algorithms include methods for determining associations between sensor data and landmarks. However, these techniques increase signicantly the complexity of the solutions. Also, they do not solve the problem of actually nding landmarks in raw sensor readings. Some examples of these algorithms include GraphSLAM [Folkesson 2004] and Sparse Extended Information Filters [Thrun. 2004].

A common technique for map representation that does not suer from data associations is occupancy grid map. In fact, occupancy grid map (OGM) [START_REF] Moravec | [END_REF], Thrun 2003] is one of the most popular environmental representation tool. It maps the environment around a vehicle as a eld of uniformly distributed binary/ternary variables indicating status of grids (occupied, free or undetected). Its applications could be dated back to 1980s, when researchers utilized sonars or lidars to illustrate environments by occupancy grids. The rst idea of employing occupancy grid map originated from the work of [Elfes 1989]. Fig. 6.1 shows an example of an occupancy grid map built by a sonar. In [Thrun 2005], occupancy grid map is regarded as the result of the mapping process and is exhaustively outlined. In addition, occupancy grid map provides an unied environmental representation for integrating dierent sensor measurements (radar, lidar, vision system).

An occupancy grid map M is dened as an evenly distributed rectangular grids within a cover area F. Each grid of the map is indexed and denoted as C i,j . The sizes of all grids are equal and set as s c . The occupancy probability of a cell C i,j is determined by corresponding sensor measurements. In order to do this eciently, we usually assume that the cells are independent, which greatly simplies the mapping process. Hence, the occupancy probability is acquired as:

P i,j (O|Y ); Y = Y 1 , ..., Y k (6.1)
where Y denotes all the possible measurements that can be associated to cell C i,j .

This chapter presents a framework to create occupancy grid maps in a dynamic environment, based on a stereo vision system and a single layer lidar. Furthermore, a fusion method to take advantages of the two sensors is presented. It has to be noted that, in previous researches [START_REF] Moravec | [END_REF],Thrun 2003], occupancy grid mapping is produced in a relatively xed environment (such as rooms, buildings) for a mobile robot. After generating a global map, it is stored for afterward usage. Under these circumstances, it is reasonable to build a global map or perform SLAM, since the robot moves in a closed loop. Whereas the situation changes in the applications of intelligent vehicle. In our application, the intelligent vehicle has to self-drive successfully in a dynamic unknown urban area. Therefore, the vehicle faces an open, unknown and dynamic environment. It emphasizes the ability of mapping in dynamic environments in real time without prior information. In our works, we do not address SLAM but try to separate mapping and localization as two tasks. Hence, the occupancy grid mapping techniques that will be presented in this chapter are about building a current local map.

The chapter is organized as follows: Sec. 6.1 proposes an occupancy grid mapping method by stereo vision. In Sec. 6.2, we describe lidar based occupancy grid mapping method. At last, a fusion of lidar and stereo vision system to build occupancy grid map is presented in Sec. 6.3 In literature, range sensors, such as lidar and radar are usually used for creating occupancy grid maps. The characteristic of measuring distance directly makes occupancy grid mapping easy to be performed. Usually, under a given sensor measurement model (such as inverse sensor model [Thrun 2003]), probabilistic occupancy grid mapping is able to be quickly calculated with the measurements. While in our application, we also use stereo vision system to observe the surroundings, which oers abundant appearance information. Nevertheless, the dierent measurement characteristics of stereo vision system lead to a dierent processing approach to get occupancy grid map.

Several existing approaches are listed as follows. In [START_REF] Murray | [END_REF]], the authors regard stereo sensor as a range nder, taking the rst encountered object in each column as an occupied grid. [START_REF] Braillon | [END_REF]] rstly estimates ground plane in front of a stereo-camera, then clusters the detected points above the plane as occupied grids. Three dierent types of occupancy grids are analyzed and compared at length in [START_REF] Badino | [END_REF]], which furthermore proposes three kinds of occupancy likelihood functions modeled by Gaussian distribution. Quite similar to [START_REF] Badino | [END_REF]], the method proposed in [START_REF] Nguyen | [END_REF]] introduces an inverse sensor model for stereocamera. In [START_REF] Danescu | [END_REF]], occupancy grid map is generated from a digital elevation map after ltering out road and isle cells according to height information. In [START_REF] Perrollaz | [END_REF]], the authors directly calculate occupancy grids by several eective probabilistic occupancy grid models in obstacle u-disparity image. In addition, this method requires a pre-performed road-obstacle separation.

However, many of the aforementioned papers do not treat the problem of moving objects in dynamic environments. Relied on the motion analysis of a moving intelligent vehicle (ego-motion estimation, moving object detection and recognition) described in Chapter. 3, 4, we propose a framework of stereo vision based dynamic occupancy grid mapping in urban environments. Dynamic occupancy grid map models real environments by evenly distributed rectangular grids, which contain both occupancy and motion information. The proposed framework mainly comprises two components (motion analysis for the vehicle-itself and independent moving objects, dynamic occupancy grid mapping) within two parallel process (sparse feature points processing between two consecutive stereo image pairs, dense stereo processing). Fig. 6.2 visualizes the whole framework. For every incoming stereo image pairs, sparse image feature points are extracted and tracked in a circular manner between the current and previous image pairs. The successfully tracked feature points are used to estimate ego-motion of the vehicle itself and independent motions of surrounding moving objects. Meanwhile, dense stereo disparity is calculated from each stereo image pair. A pixel-wise moving objects segmentation is performed in U-disparity map. Finally, the dense stereo information, together with the moving information, are used to create a dynamic occupancy grid map. Since the motion analysis has been already described in previous chapters, in this section, we will emphasize the process of creating dynamic occupancy grid map based on stereo vision.

Foundations

The sensor used for drawing occupancy grid map is still the Bumblebee XB3 stereo vision system installed in our platform SetCar. The stereo vision model used here is slightly dierent to the model introduced in Sec. 2.1.1, as we also consider the pitch angle between the ground plane and the stereo vision system. The model is illustrated in Fig. 6.3. As usual, the stereo vision system is previously calibrated and rectied. Therefore, the left and right cameras are viewed as identical and modeled by classic pinhole model (f, c u , c v ), where f is the focal length, (c u , c v ) is the position of principal point. The ground plane is assumed to be a at plane under the stereo vision system, as drawn in Fig. 6.3. The stereoscopic coordinate system is assumed to be originated from O s , the middle point of the baseline. The x-y-z directions are illustrated as well in Fig. 6.3. Therefore, 3D position of a point (X s , Y s , Z s ) in the stereoscopic coordinate system can be triangulated from its projections (u l , v l ) and (u r , v r ) in the left and right image planes as:

  X s Y s Z s   =   (u l -c u ) • b/∆ -b/2 (v l -c v ) • b/∆ b • f /∆   (6.2)
where ∆ is the disparity. Considering the pitch angle between the stereoscopic system and the ground, the corrected coordinates in the world coordinate system are:

  X w Y w Z w   =   1 0 0 0 cosθ -sinθ 0 sinθ cosθ     X s Y s Z s   +   0 h 0   (6.3)
where θ is the pitch angle, h is the height between the stereoscopic system and the ground plane. 6.1.3 Ground Plane Analysis

In urban environment, ground plane (road, parking eld, etc.) can be regarded as a supporting plane to obstacles. Usually, an obstacle roots from ground plane and is perpendicular to the plane. To successfully mapping the obstacles in surrounding environment, analyzing the ground plane is crucial. The rst component of ground plane analysis is the pitch angle estimation.

System

In most of the existing vision-based occupancy grid mapping methods (e.g. [START_REF] Badino | [END_REF]], [START_REF] Nguyen | [END_REF]], [START_REF] Perrollaz | [END_REF]]) the stereoscopic system is assumed to be parallel to the ground. However, this assumption does not hold true in our practice, since in our platform, the stereo vision system is mounted with a pitch angle to the ground. The pitch angle is crucial for mapping the surrounding environment. The importance of estimating the pitch angle between the stereo vision system and the ground plane is illustrated in Fig. points are in the coordinate of stereo vision system O s , assume there is an obstacle perpendicular to the ground, its projection on the ground plane is like a shadow (drawn in purple ellipse) caused by the pitch angle. Therefore, to make an accurate grid map of the environment, the pitch angle has to be computed and compensated.

V-disparity based Pitch Angle Estimation: As we mentioned in Sec. 4.1.2.1, one attribute of V-disparity is that, if the ground is a planar plane, it would be projected as a line, as drawn in Fig. 6.5. Furthermore, the pitch angle could be computed directly from V-disparity map. Here, we show the derivation of this property. According to the model in Fig. 6.3 and Eq. 6.3, given a 3D point in the world coordinate system O w with homogeneous coordinates P = (X, Y, Z, 1) T and its image coordinates (u, v, 1) T in the left image, we have:

λ   u v 1   = M l     X Y Z 1     (6.4)
where λ is a scale factor, M l is the transformation matrix formed by 3D rigid transformation and intrinsic projection matrix:

M l =   f c u sin θ c u cos θ b/2 0 f cos θ + c v sin θ -f sin θ + c v cos θ 0 0 sin θ cos θ 0   (6.5)
From Eq. 6.4 and 6.5, we can have: 6.6) and

v = c v + f (Y + h) cos θ -Z sin θ (Y + h) sin θ + Z cos θ (
∆ = |u l -u r | = f b (Y + h) sin θ + Z cos θ (6.7)
Assume the ground around the vehicle is at as drawn in Fig. 6.3, then its equation in O w is Y = 0. Combining the ground plane equation with Eq. 6.6 and Eq. 6.7, we can deduce the following linear equation with respect to the left image:

h b ∆ = f sin θ + (v -c v ) cos θ
(6.8) Equation 6.8 shows that a horizontal ground plane in O w will be projected as a straight line in V-disparity image I(∆, v). Considering the intercept v c in v-axis when ∆ = 0, we have: Therefore, for a planar ground, the pitch angle can be deduced as:

θ = arctan( c v -v c f ) (6.10)
The estimation of pitch angle is illustrated in Fig. 6.5. It requires to nd the line generated by the ground plane at rst, and then, compute the intercept of this line to v-axis. Then, following Eq. 6.10 gets the pitch angle estimation.

In real environments, the ground plane is not ideally at and the property of "line like" projection is hard to maintain. Hence, we use a poly-line model which consists of two-line segments to approximate real ground projection. In this polyline model, two lines (l 1 , l 2 ) correspond to two segments of ground plane according to its distance to the stereoscopic system. For example, we use one line to represent the ground within 10 meters in front of the stereoscopic system, and the other line is used to model the ground plane more than 10 meters away. We use Hough transform to extract the poly-line in V-disparity image. An example is shown in Fig. 6.6. The pitch angles of the two section of ground, θ 1 , θ 2 are estimated separately according to Eq. 6.10. Sequentially estimate pitch angles: Using two-line segments is able to estimate the pitch angles of ground sections from one stereo image pair. While in our application, the occupancy grid map is sequentially computed when the vehicle is driving in real urban area. Many factors (e.g. vibration of vehicle, ascent or downhill of the road) would contribute to the sudden changes of pitch angles. We use Kalman lter [Kalman 1960] to lter out the noise and converge toward to the true values.

In our case, we assume the pitch angles θ 1 and θ 2 as constant variables polluted by Gaussian noise and a state space model is used. Taking θ 1 for example, its state transition model is: 6.11) where θ k,1 is the current pitch angle θ 1 at time k, θ k-1,1 the state at the last time, w k,1 is a zero mean Gaussian noise with variance σ w k,1 . The measurement model is: 6.12) where z k,1 is the measurement value and v k,1 is a zero mean Gaussian noise with variance σ v k,1 . Kalman lter [Kalman 1960] is able to optimally estimate variables' posteriors in a linear dynamical system and all the error terms and measurements have a Gaussian distribution. It works in a two-step process. In the prediction step, the Kalman lter predicts estimates of the current state, along with their uncertainties. When the next measurement is observed, these predictions are updated by weighted averaging in the update step, as:

θ k,1 = θ k-1,1 + w k,1 (
z k,1 = θ k,1 + v k,1 (
θk = θk -+ W (z k -θk -) (6.13)
where θk is the posterior estimation of θ k , θk -is the prediction from state transition model. z k -θk -is the dierence between measured and predicted values, W is a weight factor. The essence of Kalman lter is a framework of recursively estimating W from the uncertainties of all the errors. The principle is to give more weight to estimates with higher certainty. More details about Kalman lter can be found in [START_REF] Welch | [END_REF]].

An example of pitch angle ltering is shown in Fig. 6.7, with an experimental result for a video sequence with almost 700 frames. The blue line is the original estimation of pitch angle θ 1 , red line is the result after Kalman ltering. The initial variances of w 0,1 and v 0,1 are set as 0.0005 and 0.001, respectively. The Kalman lter conguration for θ 2 is the same as for θ 1 . It is clear that after Kalman ltering, the pitch angle estimation is more smooth than before. Ground plane is a free space in an environment. To better map obstacles, ground areas are detected and subtracted before mapping. Detecting and subtracting ground areas are based on the tted poly-line in V-disparity map as in Sec. 6.1.3.1. As drawn in Fig. 6.8, the area below the poly-line represents ground eld. Given a pixel p(∆, v) in V-disparity image, we compute its v-distance d v to the estimated poly-line as

d v = v -v g
, where v g is its projection in the poly-line. The classication of ground and non-ground is judged by: 6.14) where d h is a threshold that restricts the ground region in V-disparity image. An example of the dened ground region is shown in Fig. 6.8. The area under the red poly-line in Fig. 6.8 (a) represents the ground region. Fig. 6.8 (b) shows the result of ground subtraction. In the process of computing occupancy probability, grids from the ground segments are directly labelled as free areas. In this section, we will describe how to build an occupancy grid map based on stereo vision system. An occupancy grid map M is dened as an evenly distributed rectangular grids C i,j array within a predened area F, where i and j are indices of a grid in the map. The sizes of all grids are equal and set to s c . Every grid C i,j holds a ternary occupancy indicator O i,j for three states:

p(∆, v) ∈ { ground if d v > d h non-ground if d v ⩽ d h (
O i,j =        undetected occupied f ree (6.15)
The process of building occupancy grid map is to estimate the status of the occupancy indicator from stereo measures. The probability of a cell C i,j to be occupied given stereo observations Y stereo is dened as P i,j (O|Y stereo ). In the following, an approach of computing P i,j (O|Y stereo ) is given. The benet of pitch angle correction is illustrated in Fig. 6.9 (we show the ground for further demonstrative purpose). Fig. 6.9 (a) -(b) show a 3D reconstruction example from a stereo image pair. The projection of all the 3D points in the ground plane without and with pitch angle correction is shown in Fig. 6.9 (c) (d), respectively. The intensity of a pixel in (c) or (d) denotes the number of points in one grid. If a grid has more points, it has higher gray value. It is clear to see that the regions belonging to the obstacles (in yellow box) become much more "brighter" after correction. This is because after correction, the reconstructed 3D points located on the surface of the obstacle become more vertical to the ground and hence more concentrated in one area. The reconstructed and corrected 3D points within ROI are assigned to corresponding grids in map M with respect to their positions. In each grid C i,j , the number of assigned points n i,j is counted.

Assuming that all obstacles are perpendicular to the planar ground, the more the number of points a grid holds, the bigger the probability it is occupied. Conse-quently, the occupancy probabilities P i,j (O|Y stereo ) are closely related to the number of associated 3D points in each grid.

However, the distribution of the associated number of 3D points is not uniform in the eld of view: the places more close to the stereoscopic system has more stereo measurements. As shown in Fig. 6.9, the ground near the stereoscopic system usually contains more associated 3D points than a real obstacle far away from the stereoscopic system. Estimating occupancy probability directly from absolute number of points would always lead to false decisions. In literature, [START_REF] Badino | [END_REF]] and [START_REF] Nguyen | [END_REF]] don't mention this problem, while [Perrollaz 2010] and [START_REF] Danescu | [END_REF]] avoid it by a previous separation of road pixels. In our method, this problem is solved by sigmoid function based adjustment described in Sec. 4.1. The absolute number of points for grid C i,j is adjusted as:

n ′ i,j = n i,j • S(d i,j ) = n i,j • r 1 + e d i,j * c (6.16)
where n i,j and n ′ i,j are the absolute and adjusted numbers of points in grid C i,j respectively. d i,j is the distance from stereo vision system to the grid. r and c are control coecients. The occupancy probability with respect to the number of points is modeled as:

P i,j (O|Y stereo ) = 1 -e -(n ′ i.j /δn) (6.17)
where δ n is a scale factor. However, it is not convenient to directly use the probability in decision making. The log-odds of the probability is then adopted: l i,j (O|Y stereo ) = log( P i,j (O|Y stereo ) 1 -P i,j (O|Y stereo ) ) (6.18) Based on the log-odds of each grid C i,j , the occupancy indicator O i,j is decided as:

O i,j =        undetected if n ′ i,j < n t occupied if l i,j (O) ≥ l t f ree if l i,j (
O) < l t or O i,j is located in the ground segments (6.19) where n t and l t are thresholds manually set for making decision. 6.1.4.3 Moving Ob ject Detection and Recognition in Occupancy Grid

Map

In chapter 4, a framework of ego-motion analysis, moving object detection, segmentation and recognition is presented. With the processing results of moving object detection and recognition, the occupancy grid map is enriched with more information. For a grid C i,j in grid map, a state vector replaces the former occupancy indicator: S i,j = (O i,j , M i,j , T i,j ), where M i,j is a binary variable indicating whether the grid is moving or not, T i,j represents the three possible types of moving objects, i.e. pedestrian, vehicle and Motor/Bike. The values of M i,j and T i,j are mainly decided by the results in chapter 4.

To label the independent grids in occupancy grid map, we distinguish the reconstructed 3D points from static environment and segmented moving objects. The 3D points from static scenes are accumulated in the grid as before and their number is represented by n s i,j . Similarly, n d i,j is the number of points from dynamic objects, we have:

n i,j = n s i,j + n d i,j (6.20)
The motion indicator M i,j is decided by comparing n d i,j with n s i,j .

M i,j = { dynamic if n d i,j > n s i,j
static otherwise (6.21) The type indicator T i,j is directly set by the recognition results from Sec. 4.2. The proposed stereo vision based dynamic occupancy grid mapping is tested using the dataset acquired by our experimental vehicle SetCar introduced in Sec. 1.3. The whole method is implemented in C++ and performed in a desktop with a CPU Intel i7-3770 quad cores 3.40GHZ. The region of interest (ROI) for the grid map is set to 30m × 30m, with a maximum height of 3m. The parameters used to calculate the occupancy indicator are set as δ n = 0.2, n t = 2, l t = 7, r = 8, c = 0.02. The whole computation time including motion analysis and dynamic occupancy grid mapping is 0.5s in average for each image pair. Four sequential results are shown in Fig. 6.10, 6.11, 6.12 and 6.13. In the gures, red areas represent detected independent moving objects, and the recognized categories are labelled outside the bounding box of the moving objects. Lidar provides an ecient approach to detect distances to the nearest obstacles. According to a certain measuring model of the sensor (e.g. inverse model [START_REF] Moravec | [END_REF]], forward model [Thrun 2003], Gaussian Beam Process [Plagemann 2007]), the occupancy probability of each cell could be quickly calculated and updated. Three common sensor models for lidar are list as follows:

• Inverse sensor model: For lidar measurements Y lidar in a cell C i,j , the probability P i,j (O|Y lidar ) species the probability of occupancy of the grid cell C i,j conditioned on the measurement Y lidar . This probability constitutes an inverse sensor model, since it maps sensor measurements back to its causes.

Occupancy grid maps usually rely on such inverse models. Notice that the inverse model does not take the occupancy of neighboring cells into account.

• Forward sensor model: Forward models proposed in [Thrun 2003] tries to overcome the assumption of inverse sensor model. It is able to calculate the likelihood of the sensor measurements for each map and set of poses, in a way that it considers all inter-cell dependencies. A forward model is of the form: P (Y lidar |M), which species a probability distribution over sensor measurements Y lidar given a map M. With forward model, mapping is transformed to an optimization problem, which aims to nd the map that maximizes data likelihood.

• Gaussian beam process: To overcome the inherent sparsity of lidar measurements, [Plagemann 2007] presented a so-called Gaussian beam process, which treats the measurement modeling task as a nonparametric Bayesian regression problem and solves it using Gaussian processes. The major benet of this approach is its ability to generalize over entire range scans directly. We can learn the distributions of range measurements for whole regions of the robot's conguration space from only few recorded or simulated range scans.

In our work, we apply the basic inverse model of lidar in occupancy grid mapping.

Ray Casting

As we described in Sec. 2.1.2, raw lidar measurements consist of detected ranges and scan angles in polar coordinate system. They can be quickly transformed into Cartesian coordinate system according to Eq. 2.1. It is reasonable to assume that every cell that lies along a line from the ego position to the detected point is free (within the maximum detected range). Before achieving this, cells traversed by a laser ray have to be marked out. We use Bresenham algorithm [Bresenham 1965] to cast rays into the grids of the map. The Bresenham algorithm is an algorithm which determines in which order one should form a close approximation to a straight line between two given points. It is commonly used to draw lines on a computer screen or raster grid. It only uses integer addition, subtraction and bit shifting, all of which are very cheap operations in standard computer architectures. Hence, its implementation is very simple and fast . The start point for the Bresenham algorithm is, in our case, always the ego car position. The endpoint is the occupied cell on that line. By doing this, we can induce all free cells corresponding to occupied regions. Areas that are occluded by occupied cells are classied as unknown. An example is show in Fig. 6.14. The gray cells between the start point and the occluded cells are all the traversed grids by the laser rays. Lidar Inverse Model

As we described previously, the probability P i,j (O|Y lidar ) represents the inverse sensor model of a 2D lidar. The objective of the inverse sensor model is to describe the formation process by which measurements are generated in a physical world and how they aect the grid cells. A lidar actively emits a signal and records it's echo. The returned echo depends on a variety of properties, such as the distance to the target, the surface material and the angle between the surface normal and the beam.

A complex lidar measuring model concerning both range and angular errors is introduced in Sec. 2.1.2, for an o-line calibration problem. However, in order to accelerate the mapping process, we simplify the measuring model by only considering http://rosettacode.org/wiki/Bitmap/Bresenham's_line_algorithm the range error. Let Y k lidar = (r k , θ k ) represent a measurement of the kth laser beam. For a cell C i,j in the grid map, its distance to Y k lidar is d(C i,j , Y k lidar ). Then, the inverse sensor model is:

P i,j (O|Y k lidar ) = λ exp(- 1 2 d(C i,j , Y k lidar ) T Σ(r k ) -1 d(C i,j , Y k lidar )) (6.22)
where λ is a scale factor. The parameter Σ(r k ) is the covariance matrix decided by r k : .23) f (r k ) reects the amount of aected regions with respect to the distance measured and is usually proportional to r k . An illustrative example of the inverse model is shown in Fig. 6.15. The radius of a point reveals its inuenced areas. Along with the increase of distance, the uncertainty of measurement grows. Hence, the inuenced areas enlarge with the distance. As mentioned before, occupancy grids not only provide information about occupied space, but also about free and unknown space. Therefore, it is necessary to dene the cells that are considered to be free or unknown given the measurement Y k lidar . Cells that are in front of the measurement should have a substantially lower occupancy probability. Inversely, cells behinds the measurement should be unknown more likely. Therefore, the nally occupancy probability is as follows:

Σ(r k ) = [ f (r k ) 0 0 f (r k ) ] ( 6
P i,j (O|Y k lidar ) = { P i,j (O|Y k lidar ) in Eq.6.15 if C i,j is in front of Y k lidar max(0.5, P i,j (O|Y k lidar )) if C i,j is behind Y k lidar (6.24)
An illustrative example of the inverse model is shown in: We implement aforementioned steps and show occupancy grid maps created by our equipped single layer lidar SICK LMS221. Angular resolution of the lidar is set to one degree, which allows to emit 181 laser beams. The maximum detected range is 80m. A grid map M covers an area of 30m × 30m. The size of each cell C i,j is 20cm×20cm. f (r k ) = r k /30, λ = 1. We only show the occupancy grid map build by lidar measures within the view of the stereoscopic system. Several results are shown in Fig. 6.17 and 6.18. In the created maps, gray grids represent free areas, black cells are unknown space and white cells are occluded regions. Around the hit point of the lidar beam, inverse lidar model spread the uncertainty due to the measured distance. From the experimental results, we can nd that, lidar measures and their corresponding occupancy grid maps are well in measuring distance. However, lidar is not able to continually observe the environment (it sparsely scans only in one layer). Some obstacles are not detected in lidar based occupancy grid map. Standard data fusion methods used for navigation tasks incorporate integration of data from sets of onboard sensors of the same type, or at least of those based on similar sensing principles. Typical classes for such sensors are range nders, data from other localization systems, etc. As the goal is to achieve maximum robustness of after fusion data, it is assumed that further improvements might be obtained by direct integration across dierent kinds of sensing systems. It is expected that this will introduce new possibilities in the fusion eld. The expectation is to bring better performance through exploring capabilities to combine multiple sources in situations when some sensors fail. A straightforward candidate for this is direct fusion of range measuring devices with intensity images. As we stated before, the lidar is superior to stereoscopic system in distance measuring. In adverse, stereoscopic system provides more integral appearance descriptions of the environment. In literature, many methods are proposed to integrate data from lidar and stereoscopic system. These methods can be roughly divided into two categories:

• Object level fusion. In [Flórez 2011], moving objects are rstly detected by a multi-layer lidar. Then, the objects are veried by a stereoscopic system to improve integrity of a driving assistance system. Similar to this, the method proposed in [START_REF] Aycard | [END_REF]] seperately detects moving objects by a multi-layer lidar and a stereoscopic system. The independently detected objects are merged for further tracking and classication.

• Map level fusion. In [Stepan 2005], occupancy grid maps are built independently by a monocular camera and a lidar. With considering dierent accuracy of the two sensors, the combination of them leads to more accurate 2D maps of the environment. [Moras 2011] proposed a credibilist approach used to model the sensor information to create occupancy grid map. Conicting information are handled by Dempster-Shafer theory. Another fusion model named linear opinion pools is proposed in [START_REF] Adarve | [END_REF]]. It models the condences of lidar measures and stereo data for integrating them in an occupancy grid map.

Fusing by Linear Opinion Pool

In our application, we choose the linear opinion pool proposed in [DeGroot 1974[START_REF] Adarve | [END_REF]] to integrate the lidar based and stereo based occupancy grid maps together. The original idea is to perform fusion as a weighted sum of sensor observations. The weight is estimated as a condence on every sensor's observation. By combining these condences, the opinion of a non reliable sensor is lowered by a low condence and hence, the result of the fusion process will depend on those reliable sensors which receive a higher weight.

In [START_REF] Adarve | [END_REF]], under the linear opinion pool, the posterior distribution over the occupancy of a cell C i,j gives the opinion of m sensors Y 1 ...Y m . Each sensor gives two quantities: its estimation for the occupancy of the cell P i,j (O|Y k ) and w k (C i,j ), a measure of the condence for such estimations. Under the assumption of mutually independent cells, the fusion of all sensory information is as follows:

P i,j (O|Y 1 ...Y m ) = α m ∑ k=1
w k (C i,j )P i,j (O|Y k ) (6.25) where α = [ ∑ k w k (C i,j )] -1 is a normalization factor for the weights. In our cases, the sensor measures either come from the lidar or stereoscopic system. Since we have already got occupancy grid maps from stereo vision system and lidar. Dierent to [START_REF] Adarve | [END_REF], where every measurement from stereoscopic system is assigned with a weight, we perform the fusion between the lidar measures and occupancy grid map made by stereoscopic system respectively. Hence, the fusion process becomes: Based on equation 6.26, for each sensor measure Y k (either from lidar or from stereoscopic system), we have to dene P (O|Y k ), the probability of a cell being occupied given the sensor information; and w k (C i,j ), the condence on the opinion.

Align Lidar with Stereoscopic System

The lidar measures and stereoscopic data are located in two dierent coordinate systems, R 3 stereo and R 3 lidar . Before applying a fusion of the data, the measurements from the two coordinate systems have to be correctly aligned.

The alignment is based on the extrinsic calibration results presented in chapter 5. After the extrinsic calibration, the 3D rigid transformation ϕ, ∆ from lidar to stereoscopic coordinate system is acquired. Hence, before the fusion, all the lidar measures are transformed into stereoscopic coordinate system according to Eq. 5.3. Next, we will discuss the models of condences.

Condences of Lidar and Stereoscopic System

The occupancy probability of a cell C i,j being occupied given the sensor measurement is given by Eq. 6.22 (lidar) or Eq. 6.17 (stereo). The condences of the two kinds of measurements are described as follows.

For the condence of a lidar measurement, consider the probability that an unexpected object is detected before the real hit of the laser beam. This is modeled by: where z * is the detected range, β is the value of occupancy probability at the hit position. After the hit, the condence decreases following a Gaussian function. An example is shown in Fig. 6.19. For the stereoscopic system, we model its measuring condence according to the distance. In fact, the condence models the fact that stereo vision works better at short distances, decreasing its accuracy when obstacles are far away. let d max be the maximum possible detected range of the stereoscopic system (d max = f ocallength * baselinelength), the measuring condence is modeled by: 1 -d 2 d 2 max (6.28) Therefore, based on the occupancy probability models in Eq. 6.22, 6.17 and condence models in Eq. 6.27 and Eq. 6.28, we apply the linear opinion pool (Eq. 6.26) to integrate data from lidar and stereoscopic system in an occupancy grid map.

Experimental Results

The fusion experiments are performed in the dataset acquired by our platform.

The image data and lidar data are rstly synchronized according to their storing time in our platform. Then, linear opinion pool based fusing is performed. In our experiments, the parameters for condence models are set as: β = 0.95, σ = 0.5, d max = 80m. Fig. 6.20, 6.21 and 6.22 demonstrate the fusion results. Because our lidar sparsely scans the environment in a single plane, sometimes, it would miss the objects due to sudden bump (Fig. 6.20 (a)), or road geometry (Fig. 6.22 (b)), or the black objects which would absorb the laser beams (Fig. 6.21 (a)). From the results, we can see that fusing the lidar and stereoscopic measurements makes the occupancy grid map more integral in observing the environment than the lidar based occupancy grid map. However, the fusion does not help the occupancy grid mapping results based on stereo vision system. The main reason is that, stereo vision system observes the environment much more densely than the single layer 2D lidar we used. Hence, most of the lidar detections within the view of stereoscopic system are already observed by vision system.

Conclusion and Future Works

In this chapter we described a framework of generating occupancy grid map from stereoscopic system, lidar and combination of both. Stereoscopic based occupancy grid mapping is achieved by rstly analyzing ground plane geometric characteristic (pitch angle) and pixel-wisely subtract ground plane in the image. Next, reconstructed 3D points are corrected by estimated pitch angles to improve the performance of occupancy grid map. Occupancy probability of each grid is calculated based on associated 3D points in each grid. In addition, independent moving objects and their categories are also labelled in the occupancy grid map by the processing results of chapter 4.

Then, we also create occupancy grid map based on a classic inverse sensor model of a single layer lidar. Since single layer lidar is not stable in observing the environment, we prove that fusing stereo information could improve the quality of occupancy grid map. However, for the stereo vision based occupancy grid map, fusing with 2D lidar data does not bring sucient improvements. The main reason of this problem is that, in our platform, a single layer 2D lidar is not robust enough to observe the environment.

For the future works, the current works could be improved by a Bayesian occupancy lter that is capable of accurately estimating and tracking the state of dynamic environments. Also, we intend to further develop tracking method in the occupancy grid map. To improve the fusion of stereo vision based and lidar based occupancy grid map, we hope the multi-layer lidar or 3D lidar could be helpful.

+D=FJAH %

Conclusions and Future Works

Conclusions

The problems addressed in this thesis concern several basic topics of perception systems in intelligent vehicles. Perception systems are the eyes of intelligent vehicles, which are in charge of environment analyzing and understanding. After a review of basic knowledges used in the thesis, we proposed a series of methods to analyze the environment.

In the rst, a stereo vision based ego-motion estimation method is presented in chapter 3. The stereo vision based visual odometry has been researched for decades. Hence, we just follow a classic framework shown in Fig. 3.1. Although this framework is relatively stable, in literature, it lacks engineering comparisons of dierent image feature detectors, as well as dierent kinds of feature association methods. We compare several representative feature detectors and proposed several standards to choose the suitable feature detectors. In addition, we also compare tracker based with matching based feature association approaches. From the comparison, we nd that CenSURE (or STAR) feature detector with KLT feature tracker achieve the best performance. The selected approach is tested both in an open dataset and our dataset. The real experimental results show that, compared with GPS, visual odometry still lacks accuracy, especially in long distance.

In the second, an independent moving object detection, segmentation and recognition method is proposed in chapter 4. The moving object detection is based on the results of visual odometry. The outliers discarded by ego-motion estimation are used for detecting moving objects in U-disparity image. The segmentation is also performed in U-disparity image by a ood-ll algorithm. Then, spatial features of segmented moving objects are extracted and sent to a pre-trained classier for recognizing three common categories of objects: pedestrian, vehicle and motor/cyclist. In our proposed method, kernel PCA algorithm is adopted to boost original simple spatial features by mapping the low dimensional data into high dimensional space. Several commonly used classiers, such as random forest, gradient boost trees are compared to get the best classier. The real experiments show that gradient boost trees achieve best recognition rate.

In the third, we developed a method to achieve optimal extrinsic calibration between the stereoscopic system and a 2D lidar using 3D reconstruction from stereo vision. We improve the accuracy of extrinsic calibration by considering the measuring models of the two sensors. The introduction of sensor models make the estimation (stereo estimation of 3D planes and nal estimation of the extrinsic pa-rameters) more precise. The improvement of accuracy is demonstrated by comparing our method to a popular method.

At last, occupancy grid map, a classic tool of interpreting environment, is built by stereoscopic system, lidar and fusing them together. For better mapping the environment by stereoscopic system, we estimate the pitch angle of the ground plane with regards to the stereoscopic system. The pitch angle compensation greatly improve the quality of occupancy grid map. Furthermore, the moving object detection and recognition results acquired previously are integrated to enhance the occupancy grid map. As for a complementary sensor to stereoscopic system, we also create occupancy grid map from 2D lidar measurements. Then, a linear opinion pool based fusing method is proposed to integrate the two kinds of measurements into one occupancy grid map. Real experimental results show that, the fusion improve the quality of lidar based occupancy grid map. However, the fusion does not improve the quality of stereo based occupancy grid map, due to 2D lidar's sparse data.

7.2

Future Works

In the author's point of view, many perspectives are envisaged to improve the current works in the future.

To improve the accuracy of stereo vision based visual odometry, several approaches could be tried:

• Develop more accurate feature detectors and feature tracking method.

• Improve the accuracy of current feature detector by estimating the sub-pixel positions.

• Consider fusing the visual odometry with other sensors, such as IMU.

For the independent moving object detection/segmentation/recognition, the possible improvements maybe:

• Combine the moving object detection with object recognition, which maybe solve the problem of failure in detecting moving pedestrians.

• Better segmentation method directly performed in raw image should be tried.

The current segmentation in disparity image is not quite satised.

• The extracted spatial features could boost the visual features (HoG, LBP, Haar). We hope to furthermore improve the recognition results by combining these two kinds of features.

For the works of calibration, in the future, we plan to introduce more dierent sensors, e.g. multi-layer lidar, 3D lidar. Hence, we plan to calibrate a multiple sensor system.

At last, for the occupancy grid mapping of the environment, the potential improvements are:

• Bayesian ltering framework could be implemented to guide the update of the map in sequence.

• Implement object tracking in the occupancy grid map.

• Semantic environmental interpretation could be added into the occupancy grid map. For example, we hope to know where is the road, buildings, pedestrians and vehicles in the map. The semantic information could be either from the object detectors or semantic image segmentation methods.

• To improve the fusing results of lidar and stereo vision system, ,ore stable and powerful lidar, such as multi-layer lidar, 3D lidar could replace the original single layer lidar.

After these possible improvements, an autonomous navigation system is able to achieved based on the current works. We hope the works presented in the thesis would be helpful for the following Ph.D students who research in this area.

R ésum é :

Les syst èmes de perception, qui sont à la base du concept du v éhicule intelligent, doivent r épondre à des crit ères de performance à plusieurs niveaux afin d'assurer des fonctions d'aide à la conduite et/ou de conduite autonome. Les travaux de cette th èse concernent le d éveloppement d'un syst ème de perception à base d'un capteur de vision st ér éoscopique et d'un capteur lidar pour l'analyse de sc ènes dynamiques en environnement urbain. Les travaux pr ésent és sont divis és en quatre parties. La premi ère partie pr ésente une m éthode d'odom étrie visuelle bas ée sur la st ér éovision, avec une comparaison de diff érents d étecteurs de primitives et diff érentes m éthodes d'association de ces primitives. Un couple de d étecteur et de m éthode d'association de primitives a ét é s électionn é sur la base d' évaluation de performances à base de plusieurs crit ères. Dans la deuxi ème partie, les objets en mouvement sont d étect és et segment és en utilisant les r ésultats d'odom étrie visuelle et l'image U-disparit é. Ensuite, des primitives spatiales sont extraites avec une m éthode bas ée sur la technique KPCA et des classifieurs sont enfin entrain és pour reconnaitre les objets en mouvement (pi étons, cyclistes, v éhicules). La troisi ème partie est consacr ée au calibrage extrins èque d'un capteur st ér éoscopique et d'un Lidar. La m éthode de calibrage propos ée, qui utilise une mire plane, est bas ée sur l'exploitation d'une relation g éom étrique entre les cam éras du capteur st ér éoscopique. Pour une meilleure robustesse, cette m éthode int ègre un mod èle de bruit capteur et un processus d'optimisation bas é sur la distance de Mahalanobis. La derni ère partie de cette th èse pr ésente une m éthode de construction d'une grille d'occupation dynamique en utilisant la reconstruction 3D de l'environnement, obtenue des donn ées de st ér éovision et Lidar de mani ère s épar ée puis conjointement. Pour une meilleure pr écision, l'angle entre le plan de la chauss ée et le capteur st ér éoscopique est estim é. Les r ésultats de d étection et de reconnaissance (issus des premi ère et deuxi ème parties) sont incorpor és dans la grille d'occupation pour lui associer des connaissances s émantiques. Toutes les m éthodes pr ésent ées dans cette th èse sont test ées et évalu ées avec la simulation et avec de donn ées r éelles acquises avec la plateforme exp érimentale v éhicule intelligent SetCar" du laboratoire IRTES-SET.

Mots

Abstract:

Intelligent vehicles require perception systems with high performances. The works presented in this Ph.D thesis concern several topics on cameras and lidar based perception for understanding dynamic scenes in urban environments. The works are composed of four parts. In the first part, a stereo vision based visual odometry is proposed by comparing several different approaches of image feature detection and feature points association. After a comprehensive comparison, a suitable feature detector and a feature points association approach is selected to achieve better performance of stereo visual odometry. In the second part, independent moving objects are detected and segmented by the results of visual odometry and U-disparity image. Then, spatial features are extracted by a kernel-PCA method and classifiers are trained based on these spatial features to recognize different types of common moving objects e.g. pedestrians, vehicles and cyclists. In the third part, an extrinsic calibration method between a 2D lidar and a stereoscopic system is proposed. This method solves the problem of extrinsic calibration by placing a common calibration chessboard in front of the stereoscopic system and 2D lidar, and by considering the geometric relationship between the cameras of the stereoscopic system. This calibration method integrates also sensor noise models and Mahalanobis distance optimization for more robustness. At last, dynamic occupancy grid mapping is proposed by 3D reconstruction of the environment, obtained from stereovision and Lidar data separately and then conjointly. An improved occupancy grid map is obtained by estimating the pitch angle between ground plane and the stereoscopic system. The moving object detection and recognition results (from the first and second parts) are incorporated into the occupancy grid map to augment the semantic meanings. All the proposed and developed methods are tested and evaluated with simulation and real data acquired by the experimental platform "intelligent vehicle SetCar" of IRTES-SET laboratory.

Keywords: Intelligent vehicle, Visual odometry, Moving objects detection and recognition, Cameras and lidar Extrinsic calibration, Dynamic occupancy grid mapping.

  Platform . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Structure of the manuscript . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Intelligent vehicles joined in DARPA grand challenge 2007 and Google driverless car

Figure 1 . 2 :

 12 Figure 1.2: Experimental platform: SetCar developed by IRTES-SET laboratory in UTBM

Figure 2

 2 Figure 2.1: Coordinate systems in our platform

Figure 2 . 3 :

 23 Figure 2.3: Pinhole camera geometry

•

  Figure 2.5: lens distortions

Figure 2 . 6 :

 26 Figure 2.6: Geometry of binocular stereo vision system

!Figure 2 . 7 :

 27 Figure 2.7: Calibrating a camera in Bouguet's Matlab camera calibration toolbox

Figure 2 .

 2 Figure 2.8: Extrinsic calibration results of a stereo vision system

Figure 2 .Figure 2

 22 Figure 2.10: Image undistortion and stereo rectication results

  Figure 2.12: The dierence between Ordinary Least Squares and Total Least Squares

  Figure 2.13: The dierence between ordinary least squares and RANSAC in 2D line tting, when facing outliers

Figure 2 .

 2 Figure 2.14: The dierence of a point on corner, edge and inside

Figure 2 .

 2 Figure 2.15: The eects of Non-Maximum Suppression in corner detection [Dominguez 2007]

Figure 2 .

 2 Figure 2.17: The eects of Laplacian of Gaussians

Figure 2 .

 2 Figure 2.18: Creating octaves of scale spaces to nd local extremas [Lowe 2004]

  e., with lter values -1 and 1. Several bi-level examples are shown in Fig.2.20. The circular lter in the left-most of Fig.2.20 is the most similar to Laplacian, but is hard to compute. The other listed lters from left to right are in an order of decreasing performance and computational cost.

Figure 2 .

 2 Figure2.20: CenSurE's bi-level lters[START_REF] Agrawal | [END_REF]

Figure 2 .

 2 Figure 2.21: SIFT descriptor [Lowe 2004]

Figure 2 .

 2 Figure 2.22: SURF descriptor [Bay 2008]

2. 4 . 3 . 1

 431 Kanade-Lucas-Tomasi Tracking (KLT tracker)

Figure 2 .

 2 Figure 2.24: A 3-dimensional k-d tree $

Figure 3 . 1 :

 31 Figure 3.1: Pipeline of sparse feature points based visual odometry

Figure 3 . 2 :

 32 Figure 3.2: Circularly association of detected feature points in consecutive image pairs

Figure 3 . 3 :Figure 3 . 4 :

 3334 Figure 3.3: Circular feature point association by KLT tracker (a) and feature matching (b)

Algorithm 3

 3 Stereo Visual Odometry based on RANSAC 1: repeat 2:

Figure 3 . 5 :

 35 Figure 3.5: Performances of repeatability for all the tested detectors

Figure 3 . 6 :

 36 Figure 3.6: Performances of uniformity (a), and speed (b) for all the tested detectors

Figure 3 . 7 :

 37 Figure 3.7: Survival rate performances for all the tested methods

Figure 3

 3 Figure 3.8: Accuracy performances for all the tested methods

Figure 3 .

 3 Figure 3.10: Four examples of stereo visual odometry experiments for KITTI dataset

Figure 3 .

 3 Figure 3.11: Quantitative translation error analyses

  Figure 3.12: The rst experimental results with our platform

  (a) Comparison between RTK-GPS and Stereo Visual Odometry (b) Translation error w.r.t path distance

Figure 3 .

 3 Figure 3.13: The second experimental results with our platform

Figure 3 .

 3 Figure 3.14: The third experimental results with our platform

Figure 3 .

 3 Figure 3.15: The fourth experimental results in our platform

  Dense disparity map and U-V disparity maps

Figure 4 . 1 :

 41 Figure 4.1: Dense disparity image and its UV-disparity images

  The inliers (red) and outliers (blue) classied by ego-motion estimation (b) The projection of outliers and inliers in U-disparity image

Figure 4 . 2 :

 42 Figure 4.2: Properties of U-disparity for obstacle detection

  (a) Sigmoid functions for intensity adjustment (b) U-disparity image before adjustment (c) U-disparity image after adjustment

Figure 4 . 3 :

 43 Figure 4.3: Intensity adjustment of U-disparity image

4. 1 . 3 Figure 4 . 4 :

 1344 Figure 4.4: Several moving object segmentation examples tested in our platform

  (a) An image of road scene (b) Moving object segmentation result (c) An image of road scene (d) Moving object segmentation result (e) An image of road scene (f) Moving object segmentation result

Figure 4 . 5 :

 45 Figure 4.5: More moving object segmentation examples

Figure 4 . 6 :

 46 Figure 4.6: The advantage of data analysis in high dimensional space

  (a) KPCA features for a cubic (b) KPCA features for a sphere

Figure 4 . 7 :

 47 Figure 4.7: The KPCA features for two dierent objects (Gaussian kernel) Feature Denition f 1 , f 2 , f 3 Cluster center in x,y,z directions f 4 , f 5 , f 6 Cluster variance in x,y,z directions f 7 Maximal height value in the cluster f 8 , f 9 , ..., f 27 Top 20 kernel shape descriptors

F

Figure 5 . 1 :

 51 Figure 5.1: Multi-sensor system: geometric relationships between the sensors coordinate systems

Figure 5 . 2 :

 52 Figure 5.2: Flow diagram of the proposed calibration method

Figure 5 . 3 :

 53 Figure 5.3: Middle point and optimal triangluation methods

Figure 5 . 4 :

 54 Figure 5.4: Estimated planes for dierent clouds of 3D points by total least squares

  Figure 5.5: Automatic lidar measurements extraction

  5.22) where Φ l and ∆ l are the ground truth values, Φ ′ l and ∆ ′ l are the estimated values. Since the results concerning with the right camera are almost the same, we only measure the error of the left camera. The results are shown in Fig.5.6 (a), (b). It's clear to see that our calibration method performs the best and tolerance of image noises is around 3 pixels.

Figure 5 . 6 :Figure 5 . 7 :

 5657 Figure 5.6: Rotation and translation errors w.r.t the noise level

  Figure 5.8: Rotation and translation errors w.r.t the noise level

Figure 5 . 9 :

 59 Figure 5.9: Rotation and translation errors w.r.t the number of poses

Figure 5 .

 5 Figure 5.10: Rotation and translation errors w.r.t the initial guess

  Figure 5.11: Comparison of the proposed method with Zhang's method [Zhang 2004]. In (c) and (d), the blue points (by the proposed method) are almost at the same height, while the red points (by Zhang's method [Zhang 2004]) are not.

Figure 5 .Figure 5 . 1 Figure 5 .Figure 5 .

 55155 Figure 5.12: Experimental results in outdoor environments. left column shows the images and projecting lidar measurements, right column shows the corresponding lidar data, red points denote within the view of the left camera.

Figure 6 . 2 :

 62 Figure 6.2: Flow diagram of the dynamic occupancy grid mapping

Figure 6 . 3 :

 63 Figure 6.3: Geometric model of the stereo vision system

Figure 6 . 4 :

 64 Figure 6.4: The importance of estimating pitch angle

Figure 6 . 5 :

 65 Figure 6.5: Estimating pitch angle to the ground

Figure 6 . 6 :

 66 Figure 6.6: Using two line segments to approximate ground projection in real environment. The green line represents ground projection more than 10 meters away. The blue line is the projection of ground within 10 meters.

  (a) Kalman ltering results for θ1 (b) Kalman ltering results for θ2

Figure 6

 6 Figure 6.7: Kalman ltering pitch angles

Figure 6

 6 Figure 6.8: Example of ground subtraction (d h = -10)

Figure 6

 6 Figure 6.9: Eect of pitch angle correction

Figure 6 .

 6 Figure 6.10: Experimental results of dynamic occupancy grid map: independent moving objects (red), static occluded areas (white), free areas (gray), undetected areas (black).

Figure 6 .

 6 Figure 6.11: Experimental results of dynamic occupancy grid map: independent moving objects (red), static occluded areas (white), free areas (gray), undetected areas (black).

Figure 6 .

 6 Figure 6.12: Experimental results of dynamic occupancy grid map: independent moving objects (red), static occluded areas (white), free areas (gray), undetected areas (black).

Figure 6 .

 6 Figure 6.13: Experimental results of dynamic occupancy grid map: independent moving objects (red), static occluded areas (white), free areas (gray), undetected areas (black).

Figure 6 .

 6 Figure 6.14: Ray casting using Bresenham algorithm

Lidar

Figure 6 .

 6 Figure 6.15: Inverse model of the lidar

Figure 6 .

 6 Figure 6.16: Prole of lidar inverse model for a beam hit on a obstacle at 50 meters

Figure 6 .

 6 Figure 6.17: (left) road scenarios and projection of lidar measures; (right) corresponding occupancy grid maps

  P i,j (O|Y 1 ...Y m ) = α m ∑ k=1 w k (C i,j )P i,j (O|Y 1 lidar ...Y m-1 lidar , Y m stereo )(6.26)

f

  orz ∈ (z * , z max ](6.27)

Figure 6 .

 6 Figure 6.19: The condence of lidar measurement (z * = 20, β = 0.9, σ = 1)

Figure 6 .

 6 Figure 6.20: Experimental results: (top) scenarios in urban environment; (middle) occupancy grid map based on lidar; (bottom) occupancy grid map by fusing lidar and stereoscopic measurements.

Figure 6 .

 6 Figure 6.21: Experimental results: (top) scenarios in urban environment; (middle) occupancy grid map based on lidar; (bottom) occupancy grid map by fusing lidar and stereoscopic measurements.

Figure 6 .

 6 Figure 6.22: Experimental results: (top) scenarios in urban environment; (middle) occupancy grid map based on lidar; (bottom) occupancy grid map by fusing lidar and stereoscopic measurements.

  -cl és : V éhicule intelligent, Odom étrie visuelle, D étection et reconnaissance d'objets en mouvement, Calibrage extrins èque entre cam éras et Lidar, Grille d'occupation dynamique.

  2.4.1 Local Invariant Feature Detectors . . . . . . . . . . . . . . . . 30 2.4.2 Feature descriptors . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4.3 Associating feature points through images . . . . . . . . . . . 43

		.	46
	2.1	Sensor Models

2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Warp the left image according to H l and the right image according to H r

	Stereo Rectification
	Epipolar line	Epipolar line
	Base line
	Figure 2.9: Stereo rectication eects
	4.

2. Compute the fundamental matrix 3. Compute projection transformations H l and H r that map the epipoles e l and e r to innity.

  and white denote weights in the pixels, grey regions have weight equal zero. Outputs of box lters can be computed very fast by integral images[Viola 2001] in a constant time. Let D xx , D xy and D yy be the results of the box lters drawn in Fig.2.19, the determinant of Hessian matrix is approximated by:

2.19

). In the gure, (a) Gaussian second order partial derivative in y-direction (b) Gaussian second order partial derivative in xy-direction (c) Box lter approximation for partial derivative in y-direction (d) Box lter approximation for partial derivative in xy-direction Figure 2.19: Box lters as fast approximation to second order derivatives of Gaussians. [Bay 2008] black

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Circular Feature Points Detection and Association . . . . .

	+D=FJAH !
	Stereo Vision based Ego-motion
	Estimation
	Contents
	3.1

50 3.2.1 Circular Feature Point Association by KLT Tracker . . . . . 51 3.2.2 Circular Feature Point Association by Matching . . . . . . . 51 3.3 Ego-motion Computation . . . . . . . . . . . . . . . . . . . . . 54 3.3.1 3D-2D Constraint based Ego-motion Estimation . . . . . . . 54 3.4 Experiments of Stereo Visual Odometry . . . . . . . . . . . 57 3.4.1 Comparing Dierent Feature Association Approaches . . . . 57 3.4.2 Stereo Visual Odometry in Urban Environment . . . . . . . . 62 3.5 Conclusion and Future Works . . . . . . . . . . . . . . . . . . 70

  incorporate error covariance matrix of triangulated feature points into motion estimation step. For this reason, they ob-

	Image Sequence
	Feature Points
	Detection
	Establishing 2D Feature
	Point Correspondences
	(Matching or Tracking)
	Motion Estimation
	(2D-2D/3D-2D/3D-3D)
	Optional	Local Optimization
		(Bundle Adjustment)
	rigid motion parameters

  Tomasi (KLT) tracker (described in Sec.2.4.3.1) is the most popular point tracking method in the computer vision eld. KLT tracker is able to accurately track a point, especially a corner point, by directly using spatial intensity information. In circular feature point association, feature points (e.g. Harris corners) are rstly detected in the the current left image. Then, they are tracked in two dierent ways to formulate a circle. Let p t l be a feature point in the current left image. The rst way tracks p t l in the current right image to p t r , then in previous right image to p t-1 r and nally in the previous left image to p t-1 l ; the second way directly tracks p t l in the previous left image to p ′t-1

	3.2.1	Circular Feature Point Association by KLT Tracker
	Kanade-Lucas-l the displacement between p t-1 l and p ′t-1 l	from the current left image. If

  3.3 Ego-motion ComputationAfter the previous steps, sparse feature points are extracted and their correspondences are established through four images in consecutive acquiring times. Hence, for a being processed stereo image pair at time t, we get a set of associated feature points as: S t = {s t

	1 , s t 2 , ..., s t i , ..., s t N }, where the subset s t i |i = 1, ..., N

Table 4 . 2

 42

: Extracted Features

  Gradient Boosting Trees: Gradient boosting is another type of ensemble learning method. It builds the model in a stage-wise fashion like other boosting methods do, and it generalizes them by allowing optimization of an arbitrary dierentiable loss function. A pseudo code of general gradient boosting is given by Algorithm 4:

	Algorithm 4 Generic gradient boosting

1: Input: training set {(x i , y i ), i = 1, ..., n}, a dierentiable loss function L(y, F (x)), number of iterations K 2: Initialize the model with a constant value:

  An interactive labeling tool is developed to label the category of a moving object in the dataset. In all, 2191 samples from 3 categories (Vehicle, Pedestrian, Motor/Bike) are collected. Within the dataset, nearly 70% of samples are treated as training samples, almost 30% are taken as testing samples, see in Table.4.3. All the testing samples are from our dataset, while the training samples are from both KITTI dataset and our dataset.

	Training samples	Testing samples	Total
	Vehicle	696	298	994
	Pedestrian	363	156	519
	Motor/Bike	380	163	543
	Overall	1439	617 2056
		Vehicle 48%
	Pedestrian 25%		Motor/Bike 27%

Table 4 .

 4

	3: Experimental dataset

  5.1. As in Sec. 2.1.1 and Sec. 2.1.3.1, we set R 3 lef t , R 3 right , R3lidaras left camera, right camera and lidar coordinate system, R 2 lef t and R 2 right as left image, right image coordinate system, respectively. Let the left camera coordinate system be the stereo vision coordinate system: R 3

	stereo

Table 5 .

 5 1: Extrinsic Parameters and Precision Measure Calculated in Real Data

http://en.wikipedia.org/wiki/Type_I_and_type_II_errors

http://en.wikipedia.org/wiki/Hough_transform

Experimental Results

The proposed algorithm has been implemented with Matlab and tested using both computer simulations and real data sets. Moreover, a comparison between our approach and a popular camera/lidar calibration method [Zhang 2004] is given.

Computer Simulations

Ground truth of the 3D rigid transformation (Φ l , ∆ l ) is set by rotation vector (in meters). We set several virtual 3D planes and select some points from each plane, project them into the images of the left and right cameras. The 3D points whose projections are within the left and right images are collected. Virtual lidar measurements on calibration board planes are acquired at the intersection of the chessboard plane and the scan plane of the lidar.

Publications

Journal: