FIELD OF THE INVENTION

The present disclosure relates to simultaneous localization and mapping (SLAM). In particular, it relates to SLAM methods and systems with reduced computational complexity.
BACKGROUND OF THE INVENTION

Selflocalization and mapping is one of the most basic capabilities of autonomous robots, and has received increasing attention in tandem with the development of hardware, such as smarter sensors and faster processors. However, 3D Mapping algorithms are difficult to apply directly to microrobot systems, especially Micro Unmanned Aerial Vehicles (MUAV) due to the following reasons. First, different from ground robots, MUAV are unable to provide odometer information to correct unacceptable incremental drift in GPSdenied environment. Second, microrobot systems can only provide limited computational resources, such as ultralow power embedded processors, due to payload and power limitations. Third, dense maps, which are crucial for higherlevel applications, such as collisionfree motion planning, object detection and scene understanding, require more computational resources to produce. The motivation of the present disclosure is to develop a light 3D dense mapping system that can be carried out by microrobot systems, and achieves faster computational speed with sufficient accuracy.

In the simultaneous localization and mapping (SLAM) problem, data association, which is to relate the sensors' measurements with the landmarks within maps, is one of the most important and timeconsuming parts.

In existing SLAM methods, iterative solutions play an important role and have become the bottleneck.
SUMMARY OF THE INVENTION

A simultaneous localization and mapping method according to a first aspect of the present disclosure comprises: receiving a frame of three dimensional point cloud data from a three dimensional point cloud sensor coupled to a robot; receiving inertial measurement data indicative of a change in pose of the three dimensional point cloud sensor; using the inertial measurement data to estimate a rotation between the received frame of three dimensional point cloud data and a key frame of three dimensional point cloud data; applying a rotation to the received frame of three dimensional point cloud data to obtain an aligned frame of three dimensional point cloud data, the aligned frame of three dimensional point cloud data having an orientation aligned with the key frame of three dimensional point cloud data; and estimating a translation of the aligned frame of three dimensional point cloud data from the key frame of three dimensional point cloud data.

By rotating the three dimensional point cloud to align it with the key frame, the process of estimating a translation of the three dimensional point cloud data can be decoupled to spaces with lower degrees of freedom which reduces the computational complexity of the SLAM processing.

In an embodiment, the method further comprises applying an axonornetric projection to the aligned frame of three dimensional point cloud data to obtain axonornetric image frame data and wherein translation of the aligned frame of three dimensional point cloud data from the key frame of three dimensional point cloud data comprises estimating a translation of the axonornetric image frame data in an axonometric image plane.

In an embodiment the method further comprises expanding the axonornetric image frame data into a onedimensional vector and wherein estimating a translation of the axonometric image frame data in an axonornetric image plane comprises determining a cyclic shift of the onedimensional vector.

The estimation of the translation of the aligned frame of three dimensional point cloud data from the key frame of three dimensional point cloud data may be carried out by applying a correlation filter. In some embodiments, the correlation filter is a general kemelized correlation filter (GKCF). The correlation filter may also be implemented as a kernelized correlation filter (KCF).

The method may comprise creating a new key frame if a matching criterion is below a first threshold. The method may further comprise refining the key frame using the aligned frame of three dimensional point cloud data if the matching criterion is above a second threshold. In some embodiments the second threshold is greater than the first threshold. The first matching criterion and/or the second matching criterion may be a peak to sidelobe ratio. In some embodiments the method further comprises refining the key frame using the aligned frame of three dimensional point cloud data.

According to a second aspect of the present disclosure a simultaneous localization and mapping system is provided. The system comprises: a three dimensional point cloud sensor; an inertial sensor; and a processing module configured to: receive a frame of three dimensional point cloud data from the three dimensional point cloud sensor; receive inertial measurement data indicative of a change in pose of the three dimensional point cloud sensor; use the inertial measurement data to estimate a rotation between the received frame of three dimensional point cloud data and a key frame of three dimensional point cloud data; applying a rotation to the received frame of three dimensional point cloud data to obtain an aligned frame of three dimensional point cloud data, the aligned frame of three dimensional point cloud data having an orientation aligned with the key frame of three dimensional point cloud data; and estimate a translation of the aligned frame of three dimensional point cloud data from the key frame of three dimensional point cloud data.

The three dimensional point cloud sensor may be a RGBD camera, and/or a Stereo Camera and/or a 3D laser scanner.
BRIEF DESCRIPTION OF THE DRAWINGS

In the following, embodiments of the present invention will be described as nonlimiting examples with reference to the accompanying drawings in which:

FIG. 1 shows a robot device which incorporates a SLAM system according to an embodiment of the present invention;

FIG. 2 is a flow diagram illustrating the framework for localization and dense mapping in embodiments of the present invention;

FIGS. 3a and 3b show perspective and axonometric projections;

FIG. 4 shows an example of a 6 degree of freedom point cloud rectified to align with a key frame then reprojected to axonometric images;

FIG. 5 shows the expansion of axonometric images into 1D vectors to determine translation between the images in an embodiment of the present invention;

FIGS. 6a and 6b show an axonometric depth image and its corresponding wellmatched points;

FIGS. 7a to 7d show an example of a key frame being refined using its subsequent frames;

FIGS. 8a to 8d show examples of correlation output from a sequence of three dimensional point cloud data frames;

FIG. 9 illustrates the process of selecting a new key frame based on the peak to side lobe ratio (PSR) value in an embodiment of the present invention;

FIG. 10 is a table showing a realtime accuracy comparison of a method according to an embodiment of the present invention with other methods;

FIG. 11 is a table showing performance efficiency corresponding to the data shown in FIG. 10;

FIG. 12 is a table showing average runtime of tracking and mapping on an ultralow power platform;

FIG. 13 is a table showing a realtime accuracy comparison on stereo cameras; and

FIG. 14 is a table showing an average runtime comparison for the results presented in FIG. 13.
DETAILED DESCRIPTION

The present disclosure is concerned with Simultaneous Localization and Mapping (SLAM) technology for robots, especially Micro Unmanned Aerial Vehicles (MUAV). It leverages on Attitude and Heading Reference System (AHRS) and any 3D sensors that can provide 3D point clouds such as RGBD camera, Stereo Cameras and 3D laser scanner. The technology is built upon decoupling 3D point cloud data and deriving a noniterative solution, which reduces computational requirements dramatically compared with traditional methods. Therefore, the system requires very low power and low computational capability devices that can be easily incorporated into an MUAV and provides a superior realtime and reliable dense map with high accuracy. Although designed originally for MUAV, it is envisaged that the systems and methods described herein can also be applied to Augmented Reality (AR) and other kinds of robots, such as Unmanned Ground Vehicles (UGV) and autonomous cars.

FIG. 1 is a block diagram schematically showing the features of a robot involved in the processing of sensor data to generate mapping data. The robot 100 may be a Micro Unmanned Aerial Vehicles (MUAV), alternatively, it may be different type of device such as an Unmanned Ground Vehicle (UGV) or an autonomous car. As shown in FIG. 1, the robot 100 comprises a 3D point cloud sensor 110. The 3D point cloud sensor 110 may be implemented as a RGBD camera, a Stereo Camera or a 3D laser scanner. The robot 100 further comprises an Attitude and Heading Reference System (AHRS) 120. The AHRS comprises a plurality of inertial sensors 125. The inertial sensors 125 may be arranged three axes to provide attitude information for the robot, including roll, pitch and yaw.

The robot 100 further comprises a processing module 130 which receives 3D point cloud data from the 3D point cloud sensor and attitude data from the AHRS 120.

Both the 3D point cloud and the attitude data from the AHRS 120 may be received as frames at regular time intervals. Alternatively, they may also be received in a hardware triggered manner. For example, the processing module may send a signal to the inertial sensors 125 indicating that measurements are required. Upon receiving such a signal, the inertial sensors 125 initiate at the same time and send the measurements to the processing module 130. This will make the measurements synchronized and bypass any time difference error between the measurements, resulting in a more accurate pose estimation. The processing module 130 uses the methods described in more detail below to generate a map of the surroundings and determine the location of the robot 100 within the map. The processing module 130 stores mapping information 140 indicating the map of the surroundings. The mapping information 140 may also include indications of the pose of the robot at different times. The pose of the robot may include an indication of the orientation and the location of the robot. The output of the method may include the acceleration, angular rate and any other variables that can describe the dynamics of the robot.

The mapping information may be used for localization and navigation while the estimated poses can be used for control and navigation in unmanned aerial vehicles (UAV), unmanned ground vehicles (UGV), augmented reality (AR), virtual reality (VR) and other applications.

In one embodiment, the system is implemented for a microrobot on a credit card sized UpBoard® equipped with an ultralow power mobilephonelevel processor Atom x5Z8350 with Scenario Design Power 2W. Running at 1.44 GHz with 2G RAM, this platform is very difficult for most of the existing SLAM algorithms to run in realtime. An Intel RealSense® camera is used as the 3D point cloud sensor.

FIG. 2 is a flow diagram illustrating the framework for localization and dense mapping in embodiments of the present invention. In the present disclosure, the term ‘noniterative’ is used in discussions of the process, here, this term is used to emphasize that the data association is noniterative while the processes of mapping and localization are still iterative with each other.

The proposed noniterative framework 200 shown in FIG. 2 comprises several subtasks. The process can be considered as a procedure of decouplingestimationrecoupling as follows: 3D point clouds are first decoupled into several subspaces, where data can be matched independently with lower complexity; then the camera motion is obtained from recoupling the estimation in the subspaces. Based on the single keyframe training, it is possible to predict the translation of the decoupled 1D vector using a general kernelized correlation filter (GKCF).

This process enables the fast motion estimation on ultralow power processors. The 6 DoF pose estimation is recoupled by fusing the inertial measurements including acceleration and angular rate with a simple pose estimator. Finally, 3D maps are fused with nonkeyframes by a moving average, so that the missing information of keyframes is complemented with complexity O(n).

In the method shown in FIG. 2, a pose estimation G^{+ }is used in data decoupling 210. The superscript ‘+’ indicates the posterior estimation of the pose. It can be regarded as the estimation when measurements are fused into the state. Here, the current measurements from AHRS are fused with the latest visual measurements (the RGBD images). The posterior state estimation of current pose is obtained after the latest inertial measurement is obtained. The poses of keyframes are obtained in a similar way, and as such they are all the posterior estimations. The terms ‘posterior’ and ‘prior’ are used to describe the state estimation obtained respectively ‘after’ and ‘before’ a measurement is obtained. The data decoupling 210 takes point cloud data 205 captured by the 3D point cloud sensor 110. The data decoupling 210 comprises rotation correction 212, axonometric reprojection 214 and image expanding 216. In the rotation correction 212, the point cloud data P 205 is rotated so that it is aligned with a key frame of point cloud data. Then, in the axonometric reprojection 214 an axonometric projection is applied to the aligned frame of point cloud data. The image expanding is applied to the axonometric projection of the aligned 3D point cloud data. This gives a 1D image vector x.

A correlation filter g is applied to the image vector x in the correlation output step 220. This gives a correlation output g(x) between the frame image and a key frame. The correlation output g(x) is used to determine an image translation estimation 230. The image translation estimation 230 is determined from the translation of the aligned 3D point cloud data against the key frame indicated by the translation of maximum value of correlation. This gives a translation estimation {tilde over (m)}. A depth translation estimation 240 is then carried out. This gives a 3D translation t which indicates the translation in position of the sensor on the robot.

The 3D translation t is input into a pose estimator 250 which uses the output of the Attitude and Heading Reference System (AHRS) 225 to provide a new pose estimation G. This pose estimation G^{+ }is the posterior pose estimation.

In addition to determining the motion of the robot, the method shown in FIG. 2 also updates a three dimensional map of the environment by updating the key frames and adding new key frames. This process occurs following the calculation of the correlation output 220. In step 260 a peak to side lobe ratio (PSR) of the correlation output g(x) is compared with a threshold T_{K}. If the PSR is below the threshold T_{K}, this indicates that the correlation between the current frame and the key frame is low. This triggers the introduction of a new key frame. Thus the method moves to step 265 in which a filter corresponding to a new key frame is trained. When a new keyframe is trained, the correlation filter will be calculated and stored together with its associated state estimation.

If the PSR is above the threshold T_{K}, this indicates that the existing key frame has a relatively strong correlation with the current frame. However, as described in more detail below, there may be information missing from the key frame. In step 270 the PSR is compared with a second threshold T_{M}. The threshold T_{M }is greater than the threshold T_{K }and thus indicates that there is a strong correlation between the current frame and the key frame. If the PSR is greater than this threshold T_{M}, the method moves to step 275 in which the key frame is refined using the aligned data from the current frame.

In the following sections, the method shown in FIG. 2 is described in more detail.
Data Decoupling

This section introduces the process of point cloud decoupling, from 6 Degrees of

Freedom (DoF) to 1 Degree of Freedom (DoF). As shown in FIG. 2 an update loop where the first block is also dependent on the output of last block is introduced. Therefore, data decoupling will be introduced on the assumption that a good pose estimation is available from the last update. At the time of beginning, it is just assumed to be located at the origin.
6 Degrees of Freedom (DoF) to 3 Degrees of Freedom (DoF)

To reduce the computation requirements, only the key frames are selected to represent a map. The other point clouds will be registered and fused with those keyframes, so that the missing information of the keyframes can be complemented. In this section, those nonkeyframes will be decoupled according to the keyframes.

A point cloud is a set of points with pose G∈SE(3) and are defined with specified color to represent an external surface of an object or environment. To decouple the rotational and translational DoF, we utilize the information of the posterior attitude estimation from the Pose Estimator 290 which is described in more detail below. Assume the k_{th }keyframe is denoted as K_{k}(G_{k}), a point cloud P_{i}(G_{i}) can be rotated to align with the keyframe by its orientation R^{+}∈SO(3) where the superscript+denotes the posterior estimation, so that the aligned point cloud {tilde over (P)}_{i }can be obtained through:

_{i}(
t _{i})=
R _{k} ^{+}(
R _{i} ^{+})
^{−1} (
G _{i} ^{+})

where

${G}^{+}=\left[\begin{array}{cc}{R}^{+}& {t}^{+}\\ 0& 1\end{array}\right]\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{with}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{R}^{+}\in \mathrm{SO}\ue8a0\left(3\right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{and}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{t}^{+}\in {\mathbb{R}}^{3}.$

In this sense, the original 6 DoF G_{i }is reduced to 3 translational DoF t_{i}, so that from now on, we only need to estimate the translation of the two point clouds P_{i }and K_{k}.

The posterior estimates G^{+ }and R^{+ }are obtained from fusing the motion estimation from the depth camera and inertial sensors and will be presented below.
3 DoF to 2 DoF

Computation requirements can be further reduced if the 3D translation movements can be decoupled into 3 independent axes. The principle is that the geometry properties in the 3 axes must be kept. We propose to apply axonometric projection on the aligned point cloud {tilde over (P)}
_{i}(t
_{i}) to get axonometric color and depth image I
^{C}; I
^{D}∈
^{M×N}. A 3D point p=(x; y; z)
^{T}∈P on the visible scene surface is mapped to image coordinates u=(u; v)
^{T}∈Ω through the axonometric projection model π
_{a}:
^{3} ^{2 }that is defined as:

$\left[\begin{array}{c}u\\ v\end{array}\right]={\frac{1}{{r}_{k}}\ue8a0\left[\begin{array}{c}x\\ y\end{array}\right]}_{\stackrel{\_}{c}}$

where r
_{k}∈
^{+ }is the projection resolution that can be set adaptively. It is different from the concept of focal length which is fixed. The subscript {tilde over (c)} denotes that the point coordinates are expressed in the aligned camera frame. The axonometric depth image is defined as:

(u)=[z]
_{{tilde over (c)}}

which stores the point distances to the axonometric plane versus the perspective depth image that stores the distances to the principal point. The differences between perspective projection and axonometric projection are illustrated in FIGS. 3a and 3 b.

FIGS. 3a and 3b show perspective and axonometric projections. FIG. 3a shows a perspective projection 310 of a scene comprising a quadrangular frustum pyramid 305. The prospective projection 310 projects 3D points on the principal point 312. As shown in FIG. 3a , the quadrangular frustum pyramid 305 appears as a rectangle 314 vertex of the quadrangular frustum pyramid 305 coincides with the principal point 312. FIG. 3b shows an axonometric projection 320. The axonometric projection 320 projects points on the axonometric plane. As shown in

FIG. 3b , when the quadrangular frustum pyramid 305 is projected into the axonometric plane, the size ratio of the two rectangles 322 and 324 will not change. In the present disclosure, this simple property is used for data decoupling.

The perspective projection does not keep the geometric relationship since an objects becomes smaller as its distance from the view point increases, while in axonometric projection, the original size of an object is kept. This property decouples the original 3D translation into the 2D translation on the axonometric plane and 1D translation perpendicular to the plane. Therefore, instead of estimating the movements of the point cloud in 3D space, we can now estimate the 2D translation of the axonometric images followed by 1D depth translation estimation.

Before that, one thing that needs to be considered is the projection resolution r_{k }which is physically the size of covered space by each pixel. The field of view of axonometric projection is a cuboid space with the central axis coinciding with the camera optical axis. Its covered space is determined by the image size M×N and the projection resolution r_{k}. Since points outside of the cuboid will be cropped, r_{k }is needed to be set large enough to make the cuboid cover most of the points. Meanwhile, there will be a large black margin on the image, if r_{k }is set too large, since not all the pixels can be covered by the projecting points. Moreover, the minimum cuboid envelope of a point cloud may change dramatically because of the camera movements. To be robust to the changing scenes, we implement a searching procedure listed in Algorithm 1 to set resolution r_{k }adaptively.


Algorithm 1: Pseudocode for Finding Adaptive r_{k} 



Input: Sampled Point Cloud K′_{k} 

Output: Adaptive Resolution r_{k} 

1 
Axonometric Image Size S 

2 
Cuboid c_{i }∈ C ; 
/* C is the candidates set that 


are calculated from image size S */ 

3 
c_{i}. num ← 0 

4 
for point p_{j }in K′_{k }do 

5 
find minimum Cuboid c_{i }∈ C that contains p_{j} 

6 
c_{i}. num ← c_{i}. num + 1 

7 c ← arg min_{ci}{c_{i } c_{i}.num ≥ 0.8*  K′_{k }, c_{i }∈ C} 

8 r_{k }← c/S 



Since most of the boundary points of a point cloud are outliers, according to one embodiment of the present invention, algorithm 1 tries to find the cuboid envelope that covers 80% points in the cloud. To be more efficient, it traverses a point cloud by sampling one point in every 25 points. Note that in this process, we keep the central axis of the cuboid envelope coincided with the camera's optical axis. Although the sampling is not precise, we find that the approximation is acceptable and the cost is deserved, since its complexity is only O(n=25), so that the computational time of Algorithm 1 can even be ignored compared to the calculation in the following sections. To further reduce repeated computation, r_{k }is updated by r_{k}+1 only when the keyframe K_{k+1}, is created. Hence, all the clouds {tilde over (P)}_{i }acquired between K_{k }and K_{k+1}, are projected with the resolution r_{k }and registered with the keyframe K_{k}. Another advantage of adaptive resolution is that when the camera is near to the scene surfaces, more details of the scenes are preserved with smaller r_{k}. A skilled person shall understand that other suitable percentages or sample sizes can also be used. For example, the cuboid envelope may cover any percentage in the range 0% to 100% of the points in the cloud. From experience, it has been found that values in the range 50%100% are particularly effective. The sample size is also empirical, any value up to 1n is possible (n is the number of points in a point cloud), the larger the value is, the faster the algorithm is, but the result is more coarse.

Since the point clouds are generated by pinhole cameras, when they are reprojected on the axonometric plane, some black holes may appear inside the image. This is because some invisible points for pinhole cameras are also projected on the axonometric plane. However, these black holes do not affect the translation estimation and can be filled by the map fusion process which is described below.

FIG. 4 shows an example that a 6 DoF point cloud is rectified by the attitude estimation and then reprojected to axonometric images. As shown in FIG. 4, a point cloud (a) captured by Intel RealSense is rectified by AHRS shown in (b) which is a perspective image. Then it is reprojected to an axonometric color image (c) and axonometric depth image (d). Some black holes 400 can be seen in (c) and (d), since some invisible points in (b) are also projected to the axonometric image plane.
2 DoF to 1 DoF

To estimate the translation of two axonometric images, we need to find the translation of their maximum overlap which is defined as the region of interest (ROI). To this end, we further decouple the axonometric image into a 1D vector.

Let an image:

I=[I _{(0)} , I _{(1)} , . . . , I _{(M1)}]
^{T}∈
^{M×N},

where I
^{T} _{(i)}∈
^{N }is the ith row of I, its expanded 1D vector x∈
^{MN }can be defined through an expanding operator x=ψ(I) where ψ:
^{M×N} ^{MN }is a bijective function on the condition that the matrix size (M, N) is known:

Ψ([I _{(0)} , . . . , I _{(M1)}]^{T})=[I _{(0)} ^{T} , . . . , I _{(M1)} ^{T}]^{T }

Define the vector's circular shift operator S
_{m}(x):
^{M×N} ^{MN }RMN, where m∈[0, MN−1] denotes the operator that circular shifts m elements for x, then the matrix's translation with (x, y) elements is defined as S
_{(x,y)}(·):
^{M×N} ^{MN }

_{(x,y)}(
I):=Ψ
^{−1}(
_{(m)}(Ψ(
I)))

where

$x=\{\begin{array}{cc}\left(m+1\right)/N& \mathrm{if}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\left(m+1\right)/N\ue2fa\frac{M}{2}\\ \left(m+1\right)/NM& \mathrm{else}\end{array}\ue89e\text{}\ue89ey=\{\begin{array}{cc}m\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\%\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eN& \mathrm{if}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89em\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\%\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eN\ue2fa\frac{M}{2}\\ m\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\%\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eNN& \mathrm{else}\end{array}$

The definition of matrix translation above is slightly different from our intuitive understanding. For the clarity of the presentation, an intuitive explanation is illustrated in FIG. 5 where 2D translation of a matrix can be regarded as a unique 1D circular shift of its expanded vector.

FIG. 5 shows the expansion of axonometric images into 1D vectors to determine translation between the images. As shown in FIG. 5, a first image I_{1 }is expanded into a 1D vector x_{1}=ψ(I_{1}) and the 2D translation S_{(x,y)}(I_{1}) between image I_{1 }and image I_{2 }can be mapped on to a unique circular shift of 1D vector S_{(m)}(x_{1}). In the example shown in FIG. 5, N=4, M=4, the matrix after translation I_{1}=S_{(x,y)}(I_{1}) can be obtained by substituting m=5 into the above equations for x and y hence x=1, y=1. The elements outside the ROI (that is the elements outside the rectangle 500) contribute nothing for finding translation, hence the malposition of these elements is ignored.

A side effect caused by the definition above is the malposition of the elements outside of the ROI. However, those elements contribute nothing for finding the translation, so that they can just be ignored. In practice, the expanding operation ψ(·) and its inverse operator ψ^{−1}(·) do not need to be implemented, since a dense matrix or an image is always stored as its expanded vector in the computer memory. The only thing we need to do is to reinterpret memory through a pointer, hence the complexity of decoupling the translation of an axonometric image is O(1). From now on, the problem becomes finding the 1D translation of the maximum overlap of two vectors. We will not distinguish an image and its expanding vector unless necessary.
Single KeyFrame Training

A nonlinear correlation filter is applied to predict the vector translation. Correlation filters are widely used for object tracking and alignment but have not been used for data association in SLAM problem. The filter used in this disclosure is based on a General Kernelized Correlation Filter (GKCF), and the following description gives a short description and shows its natural underlying connection to Minimum Output Sum of Squared Error (MOSSE) filters and Kernelized Correlation Filters (KCF).

For the sake of simplicity, the axonometric keyframe and nonkeyframes will be denoted as column vectors z
^{C,D}, x
^{C,D}∈
^{n}, where the superscript C, D denotes the color and depth images respectively. They will be left out in the absence of a little ambiguity and will not be distinguished unless necessary. The derived results can be easily extended to matrix case.

The Convolution Theorem says that the correlation operation becomes elementwise conjugate multiplication in Fourier domain. Denote the Fast Fourier Transformation (FFT) F(·) on a vector as {circumflex over (·)}, i.e. {circumflex over (x)}=F(x) so that the correlation of two vectors g=x*h is equivalent to g={circumflex over (x)}⊚
where the operator ⊚ and the superscript * denote elementwise multiplication and complex conjugate respectively. The bottleneck of this operation is to compute the forward and inverse FFTs, hence the complexity of the entire process has an upper bound O(n log n) where n is the number of elements in the vector. Define the nonlinear feature mapping function φ:
^{n} ^{d}, d>>n, the kernel trick is to find the inner product of feature mapping without calculating the high dimension features explicitly. Define the kernel function as K:
^{n}×
^{n} ^{d }such that:

(
x _{i} , x _{j})=φ(
x _{i})
^{T}φ(
x _{j})

Given a test image x∈
^{n }and its desired correlation output g∈
^{n}, the kernel correlation is defined as:

ĝ= _{Z}(
x)⊚{circumflex over (
h)}*

where:

_{Z}(
x)=[
(
x, z _{0}), . . . ,
(
x, z _{m1})]
^{T }

and z
_{i}∈
^{n }are the samplebased vectors to imply the training objective and are generated from the training sample z. Note that the length of filter h and correlation target g is related to the number of samplebased vector z
_{i}, i.e. g, h∈
^{n }and normally n≠m. The correlation output is transformed back into the spatial domain using the inverse FFT. Generally, when a set of training samples z
^{i }and their associated training outputs g
^{i }are given, we hope to find a filter h to satisfy the objective of the above kernel definition. Generally, g
^{i }can take any shape and will be discussed later. Training is conducted in the Fourier domain to take advantage of the simple elementwise operation between the input and output. To find the filter h that maps training inputs to the desired outputs, the sum of squared error (SSE) between the correlation output and the desired output in Frequency domain is minimized:

$\underset{{\hat{h}}^{*}}{\mathrm{min}}\ue89e\sum _{i}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\uf605{\hat{\kappa}}_{{Z}^{i}}\ue8a0\left({z}^{i}\right)\odot {\hat{h}}^{*}{\hat{g}}^{i}\uf606}^{2}+\lambda \ue89e{\uf605{\hat{h}}^{*}\uf606}^{2}$

where the second term is a regularization to prevent over fitting. Note that the correlation filter MOSSE also minimizes SSE, however the difference is that MOSSE only finds a linear correlation filter, while our objective is to find a nonlinear correlation filter based on the kernel trick which is still at a reasonable computational cost. In this sense, MOSSE is the linear case of GKCF and is equivalent to the above kernel correlation when k_{z}(x)=x. Solving the above optimization problem requires setting zero to the first derivative of the filter h*,

$\frac{\partial \phantom{\rule{0.3em}{0.3ex}}}{\partial {\hat{h}}^{*}}\ue89e\left(\sum _{i}\ue89e{\uf605{\hat{\kappa}}_{{Z}^{i}}\ue89e\left({z}^{i}\right)\odot {\hat{h}}^{*}{\hat{g}}^{i}\uf606}^{2}+\lambda \ue89e{\uf605{\hat{h}}^{*}\uf606}^{2}\right)=0.$

Since all the operations in this optimization problem are performed in elementwise, the elements of
can be found independently, and a closedform solution can be obtained:

${\hat{h}}^{*}=\frac{{\sum}_{i}\ue89e{\hat{g}}^{i}\odot {\hat{\kappa}}_{{Z}^{i}}^{*}\ue8a0\left({z}^{i}\right)}{{\sum}_{i}\ue89e{\hat{\kappa}}_{{Z}^{i}}\ue8a0\left({z}^{i}\right)\odot {\hat{\kappa}}_{{Z}^{i}}^{*}\ue8a0\left({z}^{i}\right)+\lambda}$

where the operator ÷ denotes the elementwise division. Note that a previously proposed kernelized correlation filter (KCF) which looks similar to the above equation but there are some differences. First, starting from the kernel regression model, the target of KCF is a scalar, while ours is based on the kernel correlation operation whose target is a vector or matrix. Second, KCF is derived from the minimization of regression error in time domain, its solution in Fourier domain is obtained by diagonalizing the circulant matrix with discrete Fourier transform, which means that KCF can only be applied when the training samples are the circular shifts of one patch. On the contrary, the solution presented above is obtained by the minimization of correlation outputs error which is in Fourier domain naturally. This makes our solution applicable to arbitrary training samples. More specifically, we can assume that only one training sample z∈
^{n }is available, hence the summation for the number of training samples z
^{i }can be left out and a more compact formula is obtained,

${\hat{h}}^{*}=\frac{\hat{g}}{{\hat{\kappa}}_{Z}\ue8a0\left(z\right)+\lambda}$

Furthermore, when we set the predesigned data xi as the circular shifts of the single training sample z, i.e. x_{i}=S_{(i)}(Z), so that m=n, this compact solution becomes exactly the same with the solution to KCF. One underlying problem of KCF is that it requires that the kernel functions satisfy K(x, z)=K(Mx,Mz) for any permutation matrix M. However, GKCF can be applied to any kernel function. In this sense, KCF is a special case of GKCF when the following conditions are satisfied: First, only one training sample z is available; second, the sample based vectors are the circular shifts z_{i}∈S(z). Last, the kernel functions treat all dimensions of data equally.

Based on the GKCF, the pattern of translational motion of the keyframe z is learnt and stored in the filter h. In the experiments, the following radialbasis kernel is used for calculating the kernel vector {circumflex over (k)}_{Z}(x):

(
x, z _{i})=
h(∥
x−z _{i}∥
^{2})

To compute the kernel vector efficiently, we first expand the norm:

h(∥x−z _{i}∥^{2})=h(∥x∥ ^{2} +∥z _{i}∥^{2}−2x ^{T} z _{i})

Since ∥x∥^{2 }and ∥z_{i}∥^{2 }are constant, the kernel vector can be calculated as:

$\begin{array}{c}{\kappa}_{Z}\ue8a0\left(x\right)=\ue89e[{h\left({\uf603\uf603x{z}_{1}\uf604\uf604}^{2},\dots \ue89e\phantom{\rule{0.8em}{0.8ex}},h\ue8a0\left({\uf603\uf603x{z}_{n}\uf604\uf604}^{2}\right)\right]}^{T}\\ =\ue89eh\ue8a0\left({\uf603\uf603x\uf604\uf604}^{2}+{\uf603\uf603z\uf604\uf604}^{2}{2\ue8a0\left[{x}^{T}\ue89e{z}_{1},\dots \ue89e\phantom{\rule{0.8em}{0.8ex}},{x}^{T}\ue89e{z}_{n}\right]}^{T}\right)\end{array}$

From the correlation theory, since z_{i }are the circular shifts of z:

x*z=[x ^{T} z _{1} , . . . , x ^{T} z _{n}]^{T }

This is substituted into the above equation for the kernel vector:

$\begin{array}{c}{\kappa}_{Z}\ue8a0\left(x\right)=\ue89eh\ue8a0\left({\uf603\uf603x\uf604\uf604}^{2}+{\uf603\uf603z\uf604\uf604}^{2}2\xb7x*z\right)\\ =\ue89eh\ue8a0\left({\uf603\uf603x\uf604\uf604}^{2}+{\uf603\uf603z\uf604\uf604}^{2}2\xb7{1}^{}\ue89e\left(\hat{x}\odot \hat{z}\right)\right)\end{array}$

The bottleneck of the lower equation is the forward and backward FFTs, so that the kernel vector can be calculated in complexity O(n log n).

In the embodiment described above, a General Kernalized Correlation Filter (GKCF) is used to estimate the vector translation. In an embodiment, the translation estimation may be carried out as follows using a kernalized correlation filter (KCF) as follows.

The fundamental idea is that we want to learn the pattern of an expanded key axonometric image so that the translation of a new image can be predicted directly based on the assumption that there is enough overlap between the key frame and the new frame. Assume that the N by M key axonometric image is denoted by a column vector x∈
^{n }where n=N×M. Hence the training samples can be defined as X−[x
_{0}, x
_{1}, . . . , x
_{n1}]
^{T }where x
_{i }is obtained from cyclic shifted x by i pixels. Basically, X is a circulant matrix and can be defined as X=C(x), since it can be generated by its first row x
^{T}. The training objective is to find a regression function y
_{i}=f(x
_{i}) where each sample x
_{i }is with a target y
_{i}, so that when a test sample z is given, f(z) will be corresponding to a translationrelated target. Intuitively, if y
_{i}=i, the translation of the test samples can be predicted directly. However, since only the zeroshift sample x
_{0 }is concerned, the label y
_{0 }is set as 1, while all the others are 0. In this sense, all the nonzeroshifts are considered as negative samples, which makes the regression function more distinctive.

1) Linear Case: A linear regression function is defined as f(x)=w
^{T}x, w∈
^{n}. Given the training samples X, the objective is to find the coefficient vector w by minimizing the squared error over the regression function and the target.

${w}^{*}=\mathrm{arg}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\underset{w}{\mathrm{min}}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\sum _{i=1}^{n}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\left(f\ue8a0\left({x}_{i}\right){y}_{i}\right)}^{2}+\lambda \ue89e{\uf603\uf603w\uf604\uf604}^{2}$

where λ is a regularization parameter to prevent overfitting. This is a ridge regression problem and can be solved by setting its first derivative to zero. Define y=[y_{0}, y_{1}, . . . , y_{n1}]^{T}, a closedform solution can be found in complex domain:

w*=(X ^{H} X+λI)^{−1} X ^{H} y

where H denotes the conjugate transpose. For a 480×360 axonometric image, where n=172800, it requires to compute the inversion of 172800×172800 matrix (X^{T}X+λI), which is impossible to carry out in realtime. However, it is interesting to note that, different from a traditional machine learning problem, X=C(x) is a circulant matrix. This amazing property makes the solution easy to obtain due to the following lemma.

Lemma 1 (J. F. Henriques): If X is a circulant matrix, the solution can be converted into frequency domain:

${\hat{w}}^{*}=\frac{{\hat{x}}^{*}\odot \hat{y}}{{\hat{x}}^{*}\odot \hat{x}+\lambda}$

where ̂ denotes the discrete Fourier transform and the superscript operator * is the complex conjugate, ⊚ and ÷ are the elementwise multiplication and division respectively.

Except for the Fourier transform, all the operations in the above are elementwise. Therefore, the complexity is dominated by the fast Fourier transform (FFT) O(n log(n)) where n is the number of points in the axonometric image. While the complexity of the original solution is dominated by a matrix inversion, whose complexity has a lower and upper bounds, given by matrix multiplication O(n2 log(n)) and GaussJordan elimination method O(n^{3}). For a 480×360 image, where n=172800, the complexity ratio r∈[O(n^{2 }log(n))/O(n log(n)); O(n^{3})/O(n log(n))]=[172800; 2480000000], which implies that lots of running time can be saved if the problem is solved by this Lemma.

2) Nonlinear Case: Data samples may be linearly separable in a high dimension space although they are not in the original space. Suppose _∅(·) is a high dimension feature space, such that ∅:
^{n}→
^{n }where d>>n, a kernel k is the inner product of the feature mapping:


where z is a test sample. The solution w∈
^{d }is expressed as a linear combination of training data xi in the feature space:

$w=\sum _{i=0}^{n1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\alpha}_{i}\ue89e\phi \ue8a0\left({x}_{i}\right)$

The regression function becomes:

$f\ue8a0\left(z\right)={w}^{T}\ue89e\phi \ue8a0\left(z\right)=\sum _{i=0}^{n1}\ue89e{\alpha}_{i}\ue89e\kappa \ue8a0\left({x}_{i},z\right)$

Then minimizing the original objective function is equivalent to finding the combination coefficient α=[α_{0}, α_{1}, . . . α_{n1}]^{T}. Given the training data X, the solution becomes:

α=(K+λI)^{−1} y

where K is the kernel matrix with each element k_{ij}=k(x_{i}; x_{j}). The dimension of α_depends on the number of samples that is the length of x. Fortunately, the kernel matrix K is circulant when k(x_{i}; x_{k}) is a Gaussian kernel:

$\kappa \ue8a0\left({x}_{i},{x}_{j}\right)=\frac{1}{2\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\mathrm{\pi \sigma}}^{2}}\ue89e{e}^{\frac{{\uf603\uf603{x}_{i}{x}_{j}\uf604\uf604}^{2}}{2\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\sigma}^{2}}}$

Therefore, the solution can be calculated in frequency domain with complexity O(nlog(n)):

$\hat{\alpha}=\frac{\hat{y}}{{k}^{\hat{x}\ue89ex}+\lambda}$

where k^{xx }is the first row of the kernel matrix K=C(k^{xx}). To guarantee the robustness, all the circular shifts of a sample z are tested. Define the kernel matrix K^{zx }where each element K_{ij}=k(z_{i}; x_{j}) and z_{i }is the i^{th }row of the circulant matrix C(z), from the regression function, we have

f(z)=K ^{zx}α

where

f(z)=[f(z _{0}), f(z _{1}), . . . , f(z _{n1})]^{T }

Since K^{zx }is a circulant matrix, again we have:

K ^{xz} =C(k ^{zx})

where k^{zx }is the first row of K^{zx}. Therefore, the response of the circular shifted samples can be found in the frequency domain:

{circumflex over (f)}(z)=k^{{circumflex over (z)}x}⊚{circumflex over (α)}

There are similarities between this derivation and the GKCF derivation described above.
Translation Estimation

Based on the single keyframe training, the translational pattern of a new axonometric image can be estimated directly by regarding the new frame as a test sample.
Image Translation Estimation

To make the correlation output more distinctive, g is set as binary vector, such that the single peak is located at the first element, i.e. g_{[0]}=1, where the brackets [·] are used accessing the element with index (starting from 0). The explanation is that, after computing a kernel correlation for a test image in the Fourier domain and converting back to the spatial domain, the shift output of value 1 will correspond to the shift of the test sample. Because of the noises and occlusion, it is nearly impossible to get the exact peak value 1. Instead, the location of the maximum value can be used to indicate the translation {tilde over (m)} of the test images.

$\stackrel{~}{m}=\mathrm{arg}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\underset{i}{\mathrm{max}}\ue89e\underset{{g}_{\left(i\right)}\ue8a0\left(x\right)}{\underset{\uf613}{{\left[i\right]}_{1}^{}\ue89e\left({\hat{\kappa}}_{Z}\ue8a0\left(x\right)\odot {\hat{h}}^{*}\right)}}$

By substituting this equation into the above equations for x and y and multiplying by the adaptive resolution r_{k}, the predicted translation: [t_{[0]}, t_{[1]}]^{T }of the axonometric image can be obtained:

[t _{[0]} , t _{[1]}]^{T} =r _{k} ·[{tilde over (x)}, {tilde over (y)}] ^{T }

The complexity the above equation for hi is bounded by the calculation of the kernel vector and the inverse FFT of the correlation output ĝ(x), hence we can obtain image translation estimation: [t_{[0]}, t_{[1]}]^{T }with complexity O(n log(n)).
Depth Translation Estimation

The ROI of axonometric image x is matched with keyframe z, if x is shifted back with {tilde over (m)} elements. Inspired by this, the camera motion in the depth direction can be estimated in the following equation which averages the depth differences of the matched ROI. Since not all the pixels in the ROI are projected by valid points in a cloud, we remove these pixels inside the black holes and only take the points in the set W defined in the following equation that contains all wellmatched pixels.

${t}_{\left[2\right]}=\frac{1}{\#\ue89e\ue536}\ue89e\sum _{i\in \ue536}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left({\ue532}_{\left(\stackrel{~}{m}\right)}\ue89e\left({x}_{\left[i\right]}^{\ue523}\right){z}_{\left[i\right]}^{\ue523}\right),\text{}\ue89e\ue536=\left\{i\ue85c\rho \left({\ue532}_{\left(\stackrel{~}{m}\right)}\ue8a0\left({x}_{\left[i\right]}^{\ue522,\ue523}\right){z}_{\left[i\right]}^{\ue522,\ue523}\right)<{T}_{c,d}\right\}$

where the operator # returns the number of elements in a set and ρ(·) is a general objective function (L_{1}norm in our tests). The predefined thresholds T_{c }and T_{d }are designed to eliminate the influence of the dynamic points. This is because dynamic objects naturally cannot be matched with the shifted image S_{−{tilde over (m)}}(x) either in color or depth spaces.

FIGS. 6a and 6b show an axonometric depth image and its corresponding wellmatched points. The image in FIG. 6a is an axonometric depth image captured from real time data. The image in FIG. 6b shows the corresponding wellmatched points by pixel intensities. The higher the brightness, the more confidence the matches have.

The advantage of using the above equations is that they only require the average differences of depth image which is extremely fast to be computed and all the wellmatched points contribute to the estimation, making it robust to depth noises. Therefore, the complexity of depth translation estimation is bounded by O(n). Now, the 3D translation t_{i }of the point cloud {tilde over (P)}_{i}, relative to the keyframe K_{k }is obtained separately based on the decoupled translational motion.
KeyFrame Selection

As the camera moves, the overlap between the new point cloud and the keyframe may decrease. Hence the peak value may not be distinct enough and new keyframe is needed to ensure large overlap. Since all the new point clouds are matched with keyframes, estimated error will be accumulated when a new keyframe is created. Although loop closure is able to reduce such accumulated error, it needs the camera to revisit a location which cannot always happen. Therefore, it is needed to be very careful to create a new keyframe. There are some useful works doing this. For example, dense visual odometry (DVO) uses the logarithmic determinant of error distribution covariance matrix to determine when to create a new keyframe. Parallel Tracking and Mapping (PTAM) created several conditions for inserting new keyframes: the distance to the last keyframe is large enough; at least twenty keyframes are passed. The ORBSLAM implementation of SLAM added new conditions: imposing a minimum visual change based on percentage of tracked features. However, these conditions are added in an adhoc way, and may not be suitable for all scenarios. For example, when the camera moves very fast, the condition of passed twenty frames may prevent creating enough keyframes and will make the camera lose tracking easily. Also, when the motion is gentle, there may be too many keyframes which are not necessary. Therefore, we argue that the condition for creating new keyframes should not be relevant to the number of passed frames and the maximum overlap or the percentage of tracked features but should be able to represent the matching quality, since the overlap of two frames may not be proportional to the matching quality because of existence of dynamic objects, etc. Considering the computational cost, we propose to apply a very simple criterion, i.e. peak to sidelobe ratio (PSR) P
_{sr}:
^{n} which is a measurement of peak strength. In this criterion, the correlation output g(x) is split into the peak which is the maximum value and the sidelobe which is the rest of pixels excluding the peak:

${P}_{\mathrm{sr}}\ue8a0\left(x\right)=\frac{\mathrm{max}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{g}_{\left[i\right]}\ue8a0\left(x\right){\mu}_{s}}{{\sigma}_{s}}<{T}_{K}$

where μ_{s }and σ_{s }are the mean and standard deviation of the sidelobe. A straightforward explanation is that the peak strength of the correlation output g(x) indicates the matching confidence level. The desired output with single peak has infinite peak strength. P_{sr}(x) is regarded as a trigger to insert a new keyframe x into map when it is smaller than T_{K}. The criterion is not only able to control the minimum confidence of each matching, but also save computational time, especially when the camera is kept still or moving slowly since there is no new training data is required. Only the calculation of mean and standard deviation of the sidelobe are performed in the calculation of the criterion, hence the complexity of computing P_{sr}(x) is O(n) where n is the number of pixels in the axonometric image.
Map Refinement and Fusion

As mentioned earlier, not all the pixels of axonometric images can be reprojected by valid points, some black holes may appear. Although they have no influence on the position estimation, we hope to fill those missing information by integrating the point clouds from other frames. A very simple idea is to complement the missing information of the keyframe z by the nonkeyframes x_{k }that are already matched with the keyframe z. This is useful especially when the camera is kept still or moving slowly, since the criterion for creating new keyframe is often not satisfied, and all those new frame x_{k }can be used to refine the keyframe z. Defining the weights vector w∈Dom(z), it is proposed to refine the keyframe z^{C,D }by the moving average defined in the following equations which relate to color and depth spaces respectively.

s
_{k }←W
^{z}+
_{(−m)}(W
_{k} ^{x})+e,

←(W
^{z}⊙
+
_{(−m)}(W
_{k} ^{x}⊙
))/s
_{k},

←(W
^{z}⊙
+
_{(−m)}(W
_{k} ^{x}⊙(
−t
_{[2]})))/S
_{k},

W
^{z}←W
^{z}+
_{(−m)}(W
_{k} ^{x}),

where e is a small scalar (set as 1e
^{−7}) to prevent division by 0. The elements of the weights vector w∈
^{n }represent the weight of the corresponding pixel to be fused. Each time a frame x
_{k }is acquired, its corresponding weight vector w
_{k} ^{x }is initialized as w
_{k[i]} ^{x}←{0,1} where 1 or 0 indicates whether the pixel x
_{[i] }can be seen in the original point cloud or not. This initialization is performed together with the procedure of axonometric reprojection. When the frame x
_{k }is selected as a new keyframe, the weights vector w
_{z }is initialized by the weight vector of that frame, i.e. w
^{z}←w
_{k} ^{x}. Because of the dynamic points and outliers, some valid points in the keyframe z cannot be matched with x
_{k}. Hence, we remove the unmatched points before performing the refining. To this end, the pixels of w
_{k} ^{x}k corresponding to those points are set as 0 in as in the following equation, so that they cannot be fused into the keyframes.

w_{k[i]} ^{x}←0, if i∉W

where the set W is defined above with reference to depth translation estimation and contains all matched points indexes. To obtain higher quality map the calculation of the moving averages is performed only when P_{sr}(xk)>T_{M }where T_{M }is a threshold to control the confidence of map fusion and normally T_{M}>T_{K}. Since all the operations are elementwise, the complexity of map refinement process is O(n) where n is the number of points in frame x_{k}.

FIGS. 7a to 7d show an example of a key frame being refined using its subsequent frames. FIG. 7a shows a new key cloud. FIG. 7b shows the refined key cloud during runtime. FIGS. 7c and FIG. 7d are the same part of new and refined key point cloud in FIG. 7a and FIG. 7b respectively. Note that the keyboard on the desk is nearly perfectly complemented. The moving average operation does not result in blurring, if the image translation is estimated correctly. FIG. 7 demonstrates that the missing information of keyframes can be complemented and much detail is preserved.
Pose Estimation and Data Recoupling

It has been shown that even the lowcost, lowprecision inertial measurement units (IMUs) can significantly improve the performance on the attitude estimation for visual odometry system. To simplify the estimation process, we propose to fuse the inertial measurements from the AHRS by loosely coupled methods. As described above, the AHRS comprises sensors on three axes and may be equipped with onboard processing unit to provide reliable attitude information. Since the gyroscope and acceleration biases of the IMU are complemented in the AHRS system, the state vector is only composed of the sensor position t∈
^{3 }in the inertial world frame, its velocity v
_{w}∈
^{3 }and attitude represented by quaternion q
_{w}. This yields a 10element state vector X=[q
_{w} ^{t}, t
^{T}, v
_{w} ^{T}]
^{T }and a simple kinetic model can be established.

To approximate the estimated translational covariance on the axonometric plane, we need to further explore the correlation output Ψ^{−1}(g(x)). Intuitively, the value of each element Ψ_{[i,j]} ^{−1}(g(x)) indicates the estimated confidence of the corresponding image translation (i, j). Hence the estimated covariance can be found by calculating the covariance of weighted relative translation to the peak. However, its complexity is O(n). To reduce the computational requirements, we assume that the normalized correlation output Ψ^{−1}(g (x))/Σg(x) is a 2D Gaussian mask with the center on the peak. Since the peak value of a 2D Gaussian function is 1/(2πσ_{i}σ_{j}) where σ_{i}, σ_{j }are the standard deviation, the estimated covariance σ_{p[0]} ^{2}, σ_{p[1]} ^{2 }in the axonometric plane can be approximated with complexity O(1):

${\sigma}_{p\ue8a0\left[0\right]}^{2}={\sigma}_{p\ue8a0\left[1\right]}^{2}=\frac{{r}_{k}^{2}\xb7\sum \phantom{\rule{0.3em}{0.3ex}}\ue89eg\ue8a0\left(x\right)\ue89e\phantom{\rule{0.3em}{0.3ex}}}{2\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\pi \xb7\mathrm{max}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{g}_{\left[i\right]}\ue8a0\left(x\right)}$

where the maximum value of g_{[i]}(x) has already been calculated and r_{k }is the projection resolution defined in Algorithm 1. Note that we do not need to calculate the summation Σg(x) in complexity O(n), since the mean value of the sidelobe μ_{s }is calculated so that Σg(x)=(n−1)μ_{s}+max g_{[i]}(x).

FIGS. 8a to 8d show examples of correlation output from a sequence. FIG. 8a is the first keyframe, FIG. 8b is the 10^{th }frame, FIG. 8c is the 43^{rd }frame and FIG. 8d is the 488^{th }frame. The correlation outputs are shown in the topright of the corresponding figures whose PSR is 161, 71 and 41 respectively. For better visualization, the outputs are normalized and shown in 3D space. We find that they can be approximated by a 2D Gaussian function with the centre on the peak. This assumption is used for fast computing the estimated covariance.

The examples presented in FIG. 8 show that the correlation output may be approximated by a Gaussian mask especially when P_{sr}(x)>20, since T_{M}>T_{K}>20, this approximation can be safely used. The estimated variance in depth translation can be obtained through direct calculation:

${\sigma}_{p\ue8a0\left[2\right]}^{2}=\underset{i\in \ue536}{\mathrm{Var}}\left({\ue532}_{\left(\stackrel{~}{m}\right)}\ue8a0\left({x}_{\left[i\right]}^{\ue523}\right){z}_{\left[i\right]}^{\ue523}\right)$

where W is a point set which is defined above. Since the estimation in the three translational subspaces are decoupled and assumed to be independent, the translation covariance u_{p }is a diagonal matrix and σ_{p}=diag(σ_{p[0]} ^{2}, σ_{p[1]} ^{2}, σ_{p[2]} ^{2}). Now based on the simple Kalman filter, we are able to recouple the movements to the original 3D space. Bounded by the procedure of keyframe training and translation estimation, the overall complexity for the whole framework is O(n log n).

As described above embodiments of the present invention consisting of a set of algorithms to decrease the computational requirements for the frontend of SLAM systems have been developed. They compose the basic components of the proposed noniterative framework.
Experimental Results

In the following section we evaluate the systems described above both quantitatively and qualitatively in terms of trajectory estimation, runtime efficiency and dense mapping in extreme conditions such as fast motion, dynamic scenes and featureless environments.

In all the experiments, the standard deviation of the Gaussian kernel is set as 0:2 pixels. The regularization term in the solutions, the matching difference threshold T_{c }and T_{d }for color and depth images are all set as 0.1. In experiments, we found these parameters are not sensitive to the test environments, since different scenarios and sensors are tested and the results are not much affected by different choices of these parameters. The PSR thresholds for creating new keyframe T_{K }and the condition for map fusion T_{M }are set as 50 and 100 respectively.

FIG. 9 illustrates the process of selecting a new key frame based on the PSR value in an embodiment of the present invention. As shown in FIG. 8, only when the PSR of the correlation output goes below T_{K}, a new keyframe will be selected. The point clouds will be fused with the keyframes when PSR is higher than T_{M}.

To test the performance of the proposed framework for microrobots on ultralow power processors and compare with the stateoftheart algorithms, we will use two different platforms. One is the credit card sized UpBoard® that is equipped with an ultralow power mobilephonelevel processor Atom x5Z8350 with Scenario Design Power 2W. Running at 1:44 GHz with 2G RAM, this platform is very difficult for most of the stateoftheart algorithms to run in realtime. Therefore, for comparison and visualization purpose, we also test our framework on a standard PC running Ubuntu 16.04 with an Intel Core i74700 CPU and 8G RAM. Limited by the payloads and power consumption, we choose Intel RealSense® camera to work with the UpBoard® and Microsoft Kinect® to work with the standard PC.

Realtime trajectory estimation using depth camera is demonstrated and compared with the stateoftheart methods, RGBD SLAM, ElasticFusion and ORBSLAM2 (RGBD version). All the reported performances are based on their open source implementation on the robot operating system (ROS) and compared in terms of the root mean squared error (RMSE), absolute mean error (MEAN) and median absolute error (MAE). Since ElasticFusion does not provide ROS interface, we implement a ROS wrapper to receive images through ROS communication protocol. Some of the above methods have a requirement to finalize the state estimation after the process is terminated, while this is not suitable for lowlatency robotic systems. To evaluate and compare their realtime performance, we only accept the instant pose estimation from all the above methods.

The experiment datasets were recorded in an indoor environment equipped with a high accurate Vicon motion capture system which serves to provide ground truth. Different traveling distances, speeds, dynamics and illumination conditions are covered. Except for the 6 DoF ground truth, the data include raw measurements from the RGBD camera at 30 Hz, inertial measurements at 100 Hz from myAHRS+1 which is a low cost high performance AHRS containing a 16bit 3axis gyroscope, a 16bit 3axis accelerometer and a 13bit 3axis magnetometer.

FIG. 10 is a table showing a realtime accuracy comparison on the Kinect dataset with PC platform. Bold texts and texts with underscore indicate the best two performances. ‘’ indicates that method cannot produce an estimation for whole trajectory.

As can be seen in FIG. 10, NISLAM achieves good performance on the PC platform. We notice that the accuracy of featurebased methods, ORBSLAM2 and RGBD SLAM2 vary a lot in different datasets, this may be due to the existence of specular reflection and featuresless regions. This phenomenon indicates that the featurebased methods are quite sensitive to feature outliers. It can be seen that ElasticFusion and ours (a system implementing the methods described herein) have similar pose estimation error for different datasets. However, the ICPbased method ElasticFusion requires powerful GPU, while our method is based on GKCF that is very fast and robust to object distortion and occlusion.

The corresponding performance of efficiency is given in FIG. 11. For RGBD SLAM2, it is noticed that the realtime pose estimation performance fluctuates much more than its offline performance, which is also given for comparison. Note that except for ORBSLAM2 (RGBD version), all the other methods are able to produce and fuse dense maps.

FIG. 11 shows that the proposed framework outperforms all the other methods in terms of running time (update rate). Note that the running time of our method reported in FIG. 11 contains the summation of both tracking and mapping; while that of the other methods only contain the tracking time. This is because the implementation of tracking thread is independent of the local mapping in other methods, while they are processed sequentially in ours with one thread, although they can also be implemented in parallel. The running time varies for different datasets according to the number of trainings. If landmarks change rapidly, our method has to train new targets and thus increases a little bit of running time. For ElasticFusion, the reported accuracy is obtained from slowing down the process of sending images 6 times to make sure that there is no dropping frames. We note the high accuracy performance of ORBSLAM2 (RGBD version) in some situations. ORBSLAM is implemented as a full SLAM system with loop closure detection, while our method is an odometry system for which drift is inevitable. Due to this reason, we found that the accuracy performance of ORBSLAM is better when visual loop closure is available frequently, and is good at tracking the environments with lots of visual features. However, we have shown that, the accuracy of such featurebased methods will drop dramatically in featureless and pseudofeature environments. While our GKCFbased method is robust to object distortion and occlusion. In this sense, the proposed noniterative framework provides a new solution for the frontend (data association) of visual SLAM system, that works well in environments for which traditional methods is not suitable. Similar to ORBSLAM and also other methods, we can integrate offtheshelf algorithms such as DBoW2 and OpenFabMap to detect visual loop and to close the loop by G20.

FIG. 12 is a table showing Average runtime (ms) of tracking and mapping on an ultralow power platform UpBoard® for RGBD image size 640×480. For ORBSLAM2, the mapping process is independent with tracking.

It has been shown that NISLAM requires the least computational resources. Experiments on the UpBoard ® also show that it is the only one that can run on such ultra low power platforms in realtime (about 30Hz update rate). The average running time of tracking and mapping are listed respectively in FIG. 12, where that of ORBSLAM2 is also given for comparison. It should be noted that ORBSLAM2 cannot produce and fuse dense maps, while our method takes tracking and dense mapping into one thread and still able to run in realtime. To the best of our knowledge, the proposed method takes dense mapping to the powerweightlimited microrobot systems for the first time.

FIG. 13 shows a realtime accuracy comparison on the stereo cameras. The corresponding efficiency performance is given in FIG. 14. Compared with RGBD cameras, stereo cameras have customized baseline and are insensitive to object materials, but generally produce more depth noises, especially in featureless regions. To show the generality and robustness, performance on stereo camera is demonstrated in this section. For comparison, the stateoftheart visual inertial odometry system using stereo cameras, OKVIS is implemented on the same platform mentioned above. Experimental data is recorded in the same Vicon room using the commercial stereo camera ZED with the inertial sensor myAHRS+, whose installation parameters are obtained from the open source multisensor calibration package Kalibr. Based on the software development kit (SDK) of ZED, we are able to obtain stereo images with size 672×376 and 3D point clouds at 30 Hz in ROS.

FIG. 13 presents the estimation error where our method shows comparable accuracy compared with OKVIS.

FIG. 14 shows average running time (ms) on PC platform using stereo cameras corresponding to the results presented in FIG. 13. It is noted that the proposed method requires much less running time than OKVIS in all the datasets. Note that OKVIS is based on nonlinear optimization and is implemented highly parallel. Hence, the reported data of OKVIS in FIG. 14 only indicates the running time of that function in each iteration, while they cannot be summed up simply.

Performance in stressful environmental conditions indicates the degree to which a system can function correctly. For a SLAM system, one of the most challenging scenes is the shaking of sensors that is caused by the robots' fast motion or motor vibration. To test the robustness of the proposed framework, this section presents the qualitative performance under the condition of fast shaking. To simplify the test procedure, a sensor module including a depth camera and an AHRS is shaken by human hands for 5 times with each lasting for 30 seconds. The number of tracking failures is counted to indicate the system robustness. For a comparison, a commercial Tangoenabled tablet computer from Google® is bound together with this sensor module. A Tangoenabled device is an Android device that integrates a wideangle camera, a depth sensing camera, an IMU sensor, GPS, accurate sensor timestamping and a software stack that enables motion tracking, depth sensing, etc. This tablet is equipped with a NVIDIA® Tegra K1 GPU with 192 CUDA cores supplied by 4GB of RAM. The number of tracking failures can be easily acquired from the screen. In this test, we found that our method can keep tracking fast shaking motion all the time, while Google's Tango enabled tablet loses tracking in such case. In this sense, our proposed framework is able to work in more challenging environments and robust to fast motion.

Whilst the foregoing description has described exemplary embodiments, it will be understood by those skilled in the art that many variations of the embodiments can be made within the scope and spirit of the present invention.