Adaptive Data Compression for Robot Perception

1 Introduction

Robot perception in complex and dynamic environments is reliant on the timely processing of low-level sensor data to provide a faithful and detailed impression of the current workspace. In recent years, 3D point clouds obtained from purpose-built laser range finders have become a particularly popular environment representation in robotics: they provide a rich source of information regarding 3D geometry. Unfortunately, the data represent only a sparse, oftentimes nonuniform sampling of the world with a manyfold redundancy in information provided per datum: for example, flat surfaces will often be sampled as densely as more complex objects. While the high degree of redundancy results in an unnecessarily inflated computational cost, the non-uniform, sparse nature of the data provides a challenge for standard perception tasks such as object detection and classification or surface reconstruction. Much is to be gained, therefore, by a similarly faithful, yet more compact, continuously queryable representation of 3D environments. The development of such a representation is the subject of this paper.

We restrict our attention to point clouds generated by a simple 3D data acquisition system commonly found in robotics: a single 2D laser scanner is pushed through the workspace and a point cloud is formed by aggregation. Based on this system, Section 3 describes the parameterisation adopted in this work. Our representation is based on a non-parametric method which simultaneously generates a continuous representation of the workspace from discrete laser samples and decimates the data set retaining only locally salient samples. At the heart of our work lies Gaussian Process (GP) regression based, in this case, with support from a finite, time varying region. This sliding-window approach to GP regression is outlined in Section 4. Inspired by the basic mechanisms of GP sparsification [Quiñonero-Candela and Rasmussen, 2005; Csató and Opper, 2002; Seeger et al., 2003], the data in the support region are chosen using the GP predictive distribution itself. We refer to this process as active sampling and describe it in Section 4.1. Active sampling explicitly retains data in regions of high complexity and therefore allows us to attain decimation factors in excess of two orders of magnitude without significant degradation in fidelity. Fig. 1 shows a typical example of the output of our system.

2 RelatedWorks

The requirement of robots to operate in large, unstructured environments provides ample incentive for research into suitable models of the workspaces encountered. A large body of work addresses the problem of 3D surface reconstruction using meshing techniques (see, for example, [Früh and Zakhor, 2003; Hoppe et al., 1993]) where every datum of a point cloud forms the vertex of a polygon. Heuristics are used to achieve the minimum number of vertices required for a suitable representation. Problems are encountered when the data arise from a non-homogeneous sampling of the workspace and/or coverage is incomplete - both are frequently the case in our problem domain. Alternative approaches address the problem by fitting geometric primitives such as planes to subsets of the raw data [Hähnel et al., 2003; Triebel et al., 2005]. This requires a strong prior on which primitives will represent the workspace well. GP models have also become a popular choice in recent years to model the terrain traversed by a robot. This is due to their ability to handle incomplete data in a principled, probabilistic fashion, Examples of such approaches include [Lang et al., 2007; Plagemann et al., 2008; Vasudevan et al., 2009]. However, these works do not address the data compression problem.

While we share the common goal of accurate 3D surface reconstruction with the literature above, this work bears closest relation to that of Gaussian Beam Processes [Plagemann et al., 2007]. The authors model laser range data on a per-scan basis using a GP model to regress range on bearing. However, while our model also regresses on range, we achieve an implicit model of the entire workspace through a sliding window approach and active data selection. This provides significant advantages above and beyond a mere extension of Gaussian Beam Processes to the 3D case.

Our active sampling strategy has been inspired by the body of work addressing GP sparsification, where active subsampling strategies have been used to select information-rich data through use of information theoretic criteria [Krause et al., 2008; Deisenroth et al., 2009; Seeger et al., 2003]. Our sparsification approach is particularly similar to [Seeger et al., 2003]. However, we exploit the time-sequential nature of laser data to form an exact and inexpensive predictive distribution for use in our decision criterion.

3 The Push-Broom 3D Laser System and its Configuration Space

The 3D laser system considered in this work consists of a simple 'push-broom' configuration: a 2D laser sensor is moved through space along an arbitrary path (see Fig. 2-left). A 3D point cloud is constructed by agglomeration of the individual scans. Our aim is to form an implicit representation of workspace surfaces by processing the data gathered. Specifically, we pose queries of points in (x, y, z) on the workspace's surface as range queries along arbitrary rays emanating from the sensor at a point along the sensor's trajectory. Let a beam from the laser be parameterised as a point q ∈ S × R where q = [θ, t]T . θ denotes the 1D angular position of the laser beam (not the whole sensor unit) and t denotes the timestamp of the laser scan. We describe q as a point in sensor configuration space Q as illustrated in Fig. 2 and refer to it as a particular configuration of the laser sensor. Sensor configuration space provides a natural domain for range regression since it is closely related to the state of the laser and not that of the rest of the robotic system. The mapping from sensor configuration space to range is necessarily a well formed function - our laser sensor associates a single range r with any point in Q such that

Consider now a mapping from Q to the Euclidean workspace W

where:

For every point q in Q, V (q) provides the six d.o.f. pose of the sensor's laser beam p ∈ S3 ×R3 at the time a measurement is taken - represented as roll, pitch, yaw and position. E(p, r) maps a six d.o.f. laser beam pose p and a scalar range measurement r to a single point x in 3D Euclidean space. This parameterisation naturally eschews the non-functional relation between elevation z and (x, y) location that is commonly found in terrain mapping formulations. By keeping each operation distinct, we also decouple robot trajectory estimation V (q), from that of the regression of the laser data G(q)1 . This permits the independent relaxation of the sensor trajectory (which, for example, could be in response to loop closure events).

To form an estimate of G(q) at any arbitrary position in Q we turn to Gaussian Process Regression. Given a set of measurements D = {(qi, ri)}Ni=1 at a query point q∗ we can obtain a predictive distribution p(r∗ |q∗ ,D). In the following section we present a brief summary of how this predictive distribution is obtained from a sliding window in sensor configuration space.

4 Sliding-WindowGaussian Process Regression

Gaussian Processes (GPs) provide for non-parametric probabilistic regression. A GP consists of a collection of jointly Gaussian distributed random variables and describes a distribution over latent functions underlying observations. It is fully specified by mean μ(q) and covariance k(q, q′ ) functions. In our application we are concerned with estimating the mapping G(q) corresponding to these latent functions. Given known ranges r from different configurations Q = {qi}Ni=1

1 Note that the regression will always be dependent on the ground-truth trajectory of the vehicle, but that no knowledge of this trajectory is required in our algorithm.

and a query point q∗ with corresponding unknown target range r∗ we can write:

In the stationary case, each element of K is derived from a suitably chosen distance metric d = ‖q - q′ ‖ between two corresponding points q and q′ in Q. We explicitly account for noise in the training observations r through an additive white noise process of strength σm along the diagonal entries of K2 . The derivation of the mean E[r∗ ] and covariance V[r∗ ] of the predictive distribution p(r∗ |q∗ ,D) for a deterministic μ(q) = 0 (as is commonly used [Vasudevan et al., 2009]) are standard and can be found, for example, in [Rasmussen and Williams, 2006]

Throughout this work we use a member of the Matérn class of covariance functions as advocated in [Stein, 1999]. We note, however, that many others, including non-stationary covariance functions [Lang et al., 2007], could be adapted and substituted in its place. The Matérn class is dependent on a shape parameter ν which regulates the smoothness of the interpolation. It equates to the more standard exponential covariance function as a special case when ν = 12 , and the squared exponential as ν →∞. As suggested in [Rasmussen and Williams, 2006], we explored several common choices of ν = { 12 , 32 , 52 ,∞} over a number of workspaces varying in complexity. We found that ν = 32 consistently produced accurate surface reconstructions for a variety of support set sizes and length scales. Although the smoother covariance functions ν = { 52 ,∞} performed well for simple workspaces, they seemed over constrained in complex scenarios, and vice versa for the rough covariance function ν = 12 . Thus

where σ2 p is the process noise, l is the length scale and d = ‖q-q′ ‖Q denotes the geodesic distance between q and q′ in Q, which for this sensor configuration is the ℓ2 norm.

In the absence of data, predictions using G(q) tend to the mean function μ(q). Thus far, this mean function has been specified a priori without accounting for the data. This leads to suboptimal behaviour particularly near the boundaries of the active window. Therefore, we modify our approach by introducing a stochastic mean function such that the GP models the residuals between a local model of the mean derived from data and the support set. This view of GP regression has been used successfully in robotics applications such as [NguyenTuong and Peters, 2010] and, in our case, provides substantial gains in terms of decimation factor. Due to space limitations the details of this approach have been omitted and the reader is referred to [Smith et al., 2010] and [Smith et al., 2011] for a detailed description.

2 In this paper we have used σ2m = 0.01m2 .

The quantities of data we consider render the application of a single monolithic GP infeasible - time complexity of a naïve implementation of GP regression is cubic in N, the size of the dataset D. Instead, for each prediction we enforce a fixed support window size n, formed from the closest (in terms of the ℓ1 norm from t to t∗) measurements that have been actively accepted by our algorithm to guarantee constant time operation. The time-sequential nature of laser data ensures that this support window slides across Q as the robot progresses along its trajectory.

We are interested in minimising the computational complexity of our algorithm for practical applications. Therefore, we now devise an active sampling strategy to intelligently determine a salient subset of D, D′ to include in our active window. The advantages of this are two-fold: we can achieve a significant compression of D, and we can re-use the inversion in Equations 6 and 7 across multiple predictions, thus decreasing computational cost. Where we must incorporate new measurements into the active window we update, rather than re-calculate the inversion. This lowers the overall worst case prediction cost to O(n2 ).

4.1 Active Sampling Using KL Divergence

Measurements are actively selected for inclusion in the support region of the GP regressor based on the information they provide with respect to the current model. Specifically, this decision is taken based on whether the current model's range prediction r∗ ∼ N(μg, σ2 g ) differs significantly from each measurement rm ∼ N(μm, σ2m) at q∗ . We use the KL divergence between the two distributions to ascertain if the average additional information required to specify rm as a result of using r∗ (instead of the true rm) is greater than a threshold κ. This threshold, naturally modulated by environment complexity, affects the rate of compression that is achieved. A sensitivity analysis of our method with a changing κ is presented in Section 6. The KL divergence between two 1D Gaussian distributions has a closed form solution:

On adoption of any new measurement into the active support set an introspection step analyses the closest (in terms of ℓ1 norm in Q) rm that had previously been marked as redundant to determine if this assessment has changed in light of the new information. On adoption of this measurement further introspection steps are applied for the next closest redundant rm until there are no further adoptions. This approach is demonstrated in Figure 3. Importantly, it is the introspection steps that allow both sides of discontinuities to be retained given the final sampling of the opposite side of the discontinuity. Figure 4 demonstrates the results of this reverse sweep, and the typical D′ that is stored and used for subsequent predictions instead of D. In the worst case, where there are truly complex sections of the workspace, our algorithm performs as well as a naïve implementation by using all D, while maintaining the ability to automatically subsample simple scenes. Note that this subsampling provide an information rich subset of D which is desirable for common applications such as data registration, where careful selection of measurements can increase accuracy and robustness [Rusinkiewicz and Levoy, 2001].

5 Experimental Setup

We now analyse system performance using both synthetic data and a real dataset collected from an outdoor urban environment. In particular, we investigate how a changing κ affects the compression rate ( |D′ | |D| %) of the original point cloud.

The real data used are part of the publically available New College data set [Smith et al., 2009]. Data were gathered from a two-wheeled Segway RMP~200 platform with two vertically aligned SICK LMS291-S14 laser sensors mounted on the side of the vehicle. Pose estimates were obtained using visual odometry from a forward facing stereoscopic camera [Newman et al., 2009]. Throughout, we use the Matérn class of covariance functions, a length scale l of 8units, a process variance σ2 p of 0.05m2 and an active window size of 200 measurements. These parameters were determined empirically to produce consistently accurate surface representations, for a given compression rate, as discussed in Section 4. The reconstruction error (in metres) is defined as the ℓ1 norm in range between a hold out set of the real measurements and corresponding predictions at the same point in Q. Predictions are made as the centre of the active window passes the hold out set in Q and are conditioned on the current active set.

6 Results

Figure 5 depicts typical subsampling as κ is varied. For aggressive decimation factors of 1,000 (corresponding to κ =

3.18 nats) to we achieve a modest mean error of 0.3m. With κ = 0.3nats errors are comparable with measurement precision of 0.015m and we need retain only one sixth of the original data.

In Figure 6 results are collated across a range of κ for the four exemplars, a) through d) (see Figure 7) and the two synthetic cases, e) and f) (see Figure 4). Intuitively, as the threshold is increased the amount of data retained decreases. This is accompanied by an increase in reconstruction error as depicted by the box plots in Figure 7. As one can discern from the images and CAD models, the relative positions of the curves correspond to the scene complexity: the more complex the scene, the greater the compression rate. Scenes b) and d) are the most complex, with noisy foliage, measurements of ceilings behind window panes and discontinuities as great as 5m, compared to that of the 1m discontinuity the person in c) presents. In all cases, these complex regions have been sampled most heavily - the outline of the person can be recognised in the Euclidean plot of c) in Figure 7.

7 Conclusions

This paper presents an overview of an adaptive 3D pointcloud compression algorithm first proposed in [Smith et al., 2010]. The approach leverages a Gaussian Process framework to provide a continuous representation of the implicit surfaces underlying 3D laser point clouds commonly encountered in robot perception. Adaptive data compression is achieved via an information theoretic selection criterion applied in a sliding window. The resulting algorithm decimates point clouds of simple workspaces by factors in excess of two orders of magnitude without significant degradation in fidelity. The computational complexity of the algorithm proposed is squared in the size if of the active window, which is constant, rather than cubic in the size of the data set. For a more detailed description of this work as well as an investigation into the use of non-stationary covariance functions that are used for the GP regression at the heart of the system, the reader is referred to [Smith et al., 2010] and [Smith et al., 2011].

8 Acknowledgements

This work was funded through an Engineering and Physical Sciences Research Council CASE award with BAE Systems, and by the European Commission under grant agreement number FP7-231888-EUROPA. We would also like to thank Winston Churchill for provision of the CAD models used throughout this paper.