Maxime Meilland

Andrew I Comport

On unifying key-frame and voxel-based dense visual SLAM at large scales

Keywords: frames, Voxel-grid RGB-D Sensor Warping Prediction Estimation Learning, Update

This paper proposes an approach to real-time dense localisation and mapping that aims at unifying two different representations commonly used to define dense models. On one hand, much research has looked at 3D dense model representations using voxel grids in 3D. On the other hand, image-based key-frame representations for dense environment mapping have been developed. Both techniques have their relative advantages and disadvantages which will be analysed in this paper. In particular each representation's space-size requirements, their effective resolution, the computation efficiency, their accuracy and robustness will be compared. This paper then proposes a new model which unifies various concepts and exhibits the main advantages of each approach within a common framework. One of the main results of the proposed approach is its ability to perform large scale reconstruction accurately at the scale of mapping a building.

3D Model

I. INTRODUCTION

The problem of dense real-time localisation and mapping within complex environments is a challenging problem for a wide range of applications ranging from robotics to augmented reality. The present work is undertaken as part of a French DGA Rapid project named Fraudo which requires dense localisation and mapping in real-time so as to allow path planning for a mobile robot to traverse uneven ground and surfaces autonomously. In this context it is important to perceive the 3D surfaces of the environment in realtime and for this dense techniques are the best approach. Amongst state-of-the-art dense SLAM approaches 3D pose tracking is usually performed in a similar manner, however, the underlying 3D representation of the dense model differs. The goal of this paper is therefore to develop an efficient, accurate and robust dense visual model for localisation and mapping that unifies various aspects of existing approaches.

More precisely, this article aims to investigate the problem of real-time dense localisation and mapping from an RGB-D camera. In the last five years many different approaches have emerged to tackle this challenging problem. Mostly, research in real-time tracking and mapping has focused on feature-based techniques [START_REF] Henry | Rgb-d mapping: Using depth cameras for dense 3d modeling of indoor environments[END_REF], [START_REF] Engelhard | Realtime 3d visual slam with a hand-held rgb-d camera[END_REF], [START_REF] Sturm | 3d pose estimation, tracking and model learning of articulated objects from dense depth video using projected texture stereo[END_REF], [START_REF] Stückler | Model learning and real-time tracking using multi-resolution surfel maps[END_REF]. Dense fusion techniques have previously been studied in the context of photometric rendering and off-line structure from motion [START_REF] Furukawa | Accurate, dense, and robust multiview stereopsis[END_REF]. Several techniques take the approach of minimising an error based on a voxelized 3D model [START_REF] Curless | A volumetric method for building complex models from range images[END_REF] (model-based). More recent approaches minimise both photometric and geometric errors directly using key-frames (image-based) in the sensor space [START_REF] Comport | Accurate quadrifocal tracking for robust 3d visual odometry[END_REF], [START_REF] Meilland | Super-resolution 3d tracking and mapping[END_REF]. The aim of this paper is to compare voxel and key-frame approaches to highlight the main differences and advantages. An improved approach will then be presented which exhibits better properties.

Although the main focus of the paper will not be on the 3D pose estimation technique, it will be advocated that direct M. Meilland and A. I. Comport are with CNRS-I3S Laboratory, University of Nice Sophia Antipolis surname@i3s.unice.fr pose estimation from photometric and depth images performs better [START_REF] Meilland | Super-resolution 3d tracking and mapping[END_REF] and subsequently this type of pose estimation will be considered common to all approaches (even if in many articles feature-based or purely geometric approaches are used). Instead this paper will focus on the underlying 3D model representation that is used to stock both geometry and photometric information of a large scale scene. Large scale scene reconstruction is relatively little studied in the literature due to the added complexity of managing large amounts of sensor information consistently. Nevertheless, several approaches have been proposed including [START_REF] Comport | Accurate quadrifocal tracking for robust 3d visual odometry[END_REF], which is based on a trajectory of key-frames, [START_REF] Meilland | A spherical robot-centered representation for urban navigation[END_REF] which proposes a graph of spherical RDB-D key-frames and, more recently, [START_REF] Roth | Moving volume kinectfusion[END_REF] and [START_REF] Whelan | Robust real-time visual odometry for dense RGB-D mapping[END_REF] who propose to extend volumetric approaches to large scale by using a moving volume.

A. Volumetric 3D Modelling

Classically, volumetric reconstructions were first proposed in the computer graphics literature [START_REF] Connolly | Cumulative generation of octree models from range data[END_REF], [START_REF] Chien | Generation of volume/surface octree from range data[END_REF]. The aim is to extend the nice regular properties of 2D image grids to 3D voxel space. These approaches therefore provide a good formalism that allows for the easy development of mathematical filters and tools for manipulating voxels. Volumetric approaches typically use a Signed Distance Function (SDF) [START_REF] Curless | A volumetric method for building complex models from range images[END_REF] to store a distance to a surface in a discretised voxel grid, representing the world geometry.

One advantage of this representation is its ease and efficiency in defining and applying mathematical operators. Another advantage is that it implicitly models self occlusions from different parts of the scene. Finally, it allows to reduce noise in structure measurements in a global 3D reference frame. Computationally, most modern graphics are cards also designed to be compatible.

Initially, volumetric approaches were used on a preacquired set of images and full bundle adjustment (BA) of all camera poses was performed to minimise drift. On the other hand, incremental approaches such as SLAM have been heavily studied in the robotics literature and more recently these approaches have been used to perform incremental realtime dense mapping [START_REF] Newcombe | Kinectfusion: Real-time dense surface mapping and tracking[END_REF], [START_REF] Roth | Moving volume kinectfusion[END_REF], [START_REF] Whelan | Robust real-time visual odometry for dense RGB-D mapping[END_REF]. The effects of incremental voxel integration are only beginning to be studied through these real-time approaches. SLAM approaches aim for real-time computation but are also subject to a much larger amount of drift over time. This is often compensated by loop closure approaches and/or reduced by performing local sliding window BA [START_REF] Mei | A constant time efficient stereo slam system[END_REF]. Unfortunately, voxel-based approaches have not been designed with incremental drift in mind and performing loop-closure on a volume is a complex task that would require redefining the underlying structure.

One main aspect of volumetric models is that, in general, they represent all space with voxels, whether the space is occupied or not. Take as a simple example the Signed Distance Function (SDF), which is used to represent occupied surfaces within a 1m 3 volume. Consider the memory requirements for a 512 × 512 × 512 × 64bits SDF, where each voxel contains a 16bit signed distance measurement, a 16bit weight and a 32bit RGBA colour value therefore requiring roughly 1GB of memory. In the literature these space requirements have been reduced by using octrees to reduce the resolution of the grid in empty space. Recently, this approach has been adapted to real-time dense mapping in [START_REF] Zeng | A memory-efficient kinectfusion using octree[END_REF] and [START_REF] Damen | Egocentric real-time workspace monitoring using an rgb-d camera[END_REF].

A disadvantage of this approach is that the knowledge about the camera positions, used to acquire the map, are usually completely discarded (i.e. it is independent of the sensor trajectory) and the volumetric approach is unable to easily capture the non-linear variation of the image resolution, which depends on a particular camera path. Whilst a co-variance can be stored for each voxel, the associated uncertainty becomes linearised in 3D space and the uncertainties of the sensor measurements are approximated. Finally volumetric approaches will consume too much memory to represent high image resolutions (super-resolution [START_REF] Meilland | Super-resolution 3d tracking and mapping[END_REF]) across all space (as opposed to locally around nodes in a graph).

B. Image-based Key-frame Modelling

Image-based key-frame models have been popular in the robotics literature, due to their locally accurate representation, their capability to handle incremental drift and readjust a more compact set of poses using classic loop-closure constraints in real-time [START_REF] Kümmerle | g2o: A General Framework for Graph Optimization[END_REF]. Early real-time key-frame approaches based on features [START_REF] Nistér | Visual odometry[END_REF], [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF], [START_REF] Konolige | Frameslam: from bundle adjustment to realtime visual mapping[END_REF]. The first realtime dense visual odometry approach was published in [START_REF] Comport | Accurate quadrifocal tracking for robust 3d visual odometry[END_REF] but key-frames corresponded to real-images and no structure or photometric information was fused. Recently, a key-frame approach has been published, which fuses simultaneously geometry and photometry onto a set of super-resolution keyframes [START_REF] Meilland | Super-resolution 3d tracking and mapping[END_REF] whereby the closest key-frame in the graph is often used to perform localisation. Subsequently key-frame approaches also have great noise reduction properties.

Another advantage of key-frame approaches is their memory efficiency when the viewing path is the same as the acquisition trajectory. Contrary to volumetric approaches, these models encode the raw sensor data (and sensor uncertainty) along with the trajectory associated with the image acquisition, whereas volumetric approaches discard this information. One further result of maintaining a reference trajectory is that it directly encodes the camera resolution from these vantage points since the resolution of the map is a function of the path taken to acquire it. The resolution of volumetric approaches aims to cover all vantage points equally. In terms of comparing memory efficiency, consider the extreme boundary case of a single floating point RGB-D key-frame of dimensions 640 × 480 × 64bits. It contains 16 bits for each depth, 16 bits for each weight, 32 bits for each RGBA colour and the parameters for the pose of the key-frame can be considered negligible. In this case the keyframe image only requires 2.34 MB to represent a local subset of viewpoints through novel-view synthesis (depending on the scene structure). Furthermore, this representation can be accessed in constant time, independently of the number of key-frames in the graph.

Whilst volumetric approaches benefit from efficient computation due to a regular (non-scattered) grid pattern, keyframe approaches are also able to perform computation efficiently by exploiting a 2D grid pattern in the image space.

As will be seen later, one disadvantage of key-frame approaches is that the current camera image does not always correspond to the same field of view as the closest keyframe and subsequently only have a partial overlap. This can degrade the localisation robustness since only a small area of the current image is registered. This can be solved by more sophisticated reference images such as spherical panoramas with depth [START_REF] Meilland | Dense visual mapping of large scale environments for real-time localisation[END_REF]. However, this model is not easily applicable for commonly used sensors with limited field of view.

C. An approach to Unifying

In this paper, a new approach will be proposed that aims to combine the advantages of volumetric and key-frame based approaches. From the outset, it is clear that key-frame approaches encode additional information about the acquisition phase. In particular, the acquisition trajectory is known and raw sensor data at set resolutions are stored. This raw data can easily be transformed into a volumetric representation to exhibit the properties of volumetric approaches. Nevertheless if we discard the acquisition data from the model then it will be lost completely. As such, the proposed model will maintain an underlying key-frame representation that can easily be transformed into a volumetric 3D model if required.

The aim of the new approach will therefore be to develop a multi-key-frame approach that can easily handle self occlusions and that can fuse multiple sources of data as volumetric approaches do. A simple approach to solve this problem is to simultaneously use multiple key-frames extracted from the graph, however, this would require warping the current image several times, increasing the computational cost. To avoid the redundant warping of the current image onto several overlapping reference frames, a new forward composition approach is proposed by developing a computationally efficient warp and blending function. This function will take multiple reference images and combine them efficiently in real-time to create a single predicted frame. As will be demonstrated, this allows more efficient registration with respect to only one synthesised reference key-frame.

Returning to the small example on space requirements, volumetric approaches account for the set of all 3D points implicitly in the structure (up to a given resolution). When compared to the same resolution as a single key-frame, the voxel model requires adding a 3rd dimension to the image size. For example, a key-frame image of size 512 2 becomes a volume of dimension 512 3 to represent the same resolution for that view. Of course, a single key-frame does not cover arbitrary viewpoints as a volumetric model does and it is necessary to determine how many key-frames are needed. Unfortunately this is highly dependant on the scene structure.

In the ideal case for key-frame models, the camera observes either a convex scene or a concave object. In that case the equivalent number of images to the voxel grid corresponds to the 6 images of a cube observing the scene or the object at the same resolution as the voxel grid (assuming an orthographic camera projection model). The key-frame Fig. 1. Unified image and model-based system. For both volumetric and key-frame based representations, tracking and mapping can be decomposed in four major blocks: warping, prediction, estimation and model update model requires 512 2 × 6 images plus 6 poses to represent the scene, while, volumetric approaches require 512 3 voxels minus any compression achieved with octrees. It is clear that if complex non-concave objects or non-convex scenes need to be mapped, key-frame based approaches will require an increasing number of reference views to account for the self occlusions of the scene. However, in the proposed approach, key-frames are to be chosen carefully to encapsulate the concave and convex elements of the scene, allowing a compact and efficient representation to be found.

II. DENSE REAL-TIME TRACKING AND MAPPING

A. Introduction

Consider a calibrated RGB-D sensor with a colour brightness function

I : Ω × R + → R + ; (p) → I(p, t) and a depth function D : Ω × R + → R + ; (p, t) → D(p, t), where Ω = [1, n] × [1, m] ⊂ R 2 , P = (p 1 , p 2 , . . . , p nm ) ⊤ ∈ R mn×2
⊂ Ω are pixel locations within the image acquired at time t, and n × m is the dimension of the sensor's images. It is convenient to consider the set of measurements in vector form such that I(P, t) ∈ R +nm×1 and D(P, t) ∈ R +nm×1 . Note that t and P may be omitted for clarity.

V = (v 1 , v 2 , . . . , v nm ) ⊤ ∈ R mn×3 is defined as the matrix of 3D vertices related to the surface according to the following point-depth back-projection:

v i = K -1 p i D(p i ), (1)
where K ∈ R 3×3 is the intrinsic camera matrix and p i are homogeneous pixels coordinates. V will be taken to be a 3D vertex function

V : Ω × R + → R 3 ; (p, t) → V(p, t).
The set S = {I, V, N, C} is defined to be a 3D textured surface. The uncertainties C ∈ R nm are computed for each vertex according to the depth sensor model of [START_REF] Khoshelham | Accuracy and resolution of kinect depth data for indoor mapping applications[END_REF]. The normals N ∈ R 3×nm are computed for each vertex using a local cross product on the image grid. Again, the surface normals and the uncertainties will be considered as functions such that

N : Ω × R + → R 3 ; (p, t) → V(p, t) and C : Ω × R + → R; (p, t) → C(p, t).
A typical localisation/reconstruction pipeline is shown in Figure 1. For both volumetric and key-frame based approaches, the incremental reconstruction pipeline contains the following blocks: frame prediction, warping, estimation and model update. Here frame-to-frame pose estimation is made between a predicted view and the live frame so that this step is the same for both approaches. The only blocks that differ are the view prediction and model update.

B. Pose estimation

Now consider the predicted reference surface S * = {I * , V * , N * , C * } to be a view prediction of the current surface S = {I, V, N, C}. The goal is to find the unknown motion parameters x ∈ R 6 between S * and S defined as:

x = 1 0 (ω, υ)dt ∈ se(3), (2)
which is the integral of a constant velocity twist which produces a pose T, where T = (R, t) ∈ SE(3). R ∈ SO(3) is a rotation matrix and t ∈ R(3) a translation vector. The pose and the twist are related via the exponential map as T = e [x]∧ with the operator [.] ∧ as:

[x] ∧ = [ω] × υ 0 0 , (3)
where [.] × the skew symmetric matrix operator.

An improved version of the bi-objective direct iterative closest point (ICP) approach of [START_REF] Tykkala | Direct iterative closest point for real-time visual odometry[END_REF] is employed, which simultaneously minimises a photometric error along with a geometric error between the surfaces S * and S such that

e(x) =   I w( TT(x), V * ) -I * (P * ) RR(x)N * ⊤ V ⊤ -Π TT(x)V * ⊤ ⊤   , (4)
where the first row of equation ( 4) is the photometric term and the second row is a point-to-plane ICP error with projective data association. The function w(•) is the inverse warping function described in Section III-B. This non-linear error is iteratively minimised using a Gauss-Newton approach such that:

x = -(J T WJ) -1 J T We, (5)
where J T = J esm N * ⊤J icp . J contains the stacked Jacobian matrices of the errors of equation ( 4). J esm is the Jacobian matrix of the photometric term computed using the efficient second order (ESM) approach [START_REF] Malis | Improving vision-based control using efficient secondorder minimization techniques[END_REF] and J icp is the standard ICP Jacobian matrix. W is a diagonal weighting matrix of dimensions 2nm × 2nm obtained by M-estimation [START_REF] Zhang | Parameter Estimation Techniques: A Tutorial with Application to Conic Fitting[END_REF], with Huber's influence function. Since the photometric error is influenced by the geometric error, the matrix W is computed such that:

W = W esm W icp 0 0 W icp . (6)
This weighting scheme efficiently handles geometric occlusions, such as dynamic moving objects. Note that the photometric errors (and subsequently the weights) are coupled with the geometry while the geometric error is independent of the photometry (only for projective light RGB-D sensors).

The pose estimate T is finally updated using a homogeneous update until convergence as T ← TT(x).

III. SYNTHESISING A NOVEL VIEW FROM THE MODEL

The main difference between volumetric and key-frame approaches appears whenever it is necessary to interact and use the 3D model data. Even within each of these approaches several errors may be minimised, primarily depending on a forward or inverse warping function. These warping functions will be elaborated here and a new warping function will be proposed for combining several key-frames together efficiently, evidently solving the problem of partial overlap and occlusion (between current and key-frame images).

A. Geometric warping

The geometric warping function P w⊤ = w (T, K, V), transforms the vertices V with the pose T and projects the transformed vertices on the image plane by

P w⊤ = MV ⊤ ,
where the perspective projection matrix M is defined such that M = KΠT. Π = [1, 0] ∈ R 3×4 projects the 4 × 4 pose matrix onto the 3 × 4 space and P w⊤ is normalised by dividing by its third co-ordinate to obtain P w⊤ .

B. Inverse warping / texture projection

Inverse warping is the equivalent of texture projection in computer graphics. It basically consists in projecting the current image/texture onto the world geometry. It is often referred as inverse warping since the transformation from the source intensities to the destination space is performed by projecting the destination vertices (reference image) to the source space (current image) and then interpolating the intensity values on a regular grid as:

I * (P * ) = I(w(T, K, V * )). (7)
As it can be seen in equation ( 7), the current intensity warping only depends on the reference vertices V * . This allows to perform computationally efficient grid interpolation of the surface and obtain a correspondence between the reference and the current intensities. One major drawback of this approach is that occlusions are not handled, since several vertices may project onto the same pixel, leading to inconsistencies in the warped image. This is, however, the fastest way to transfer image intensities between frames and its Jacobian can be pre-computed.

C. Forward warping / rasterisation

Standard rasterisation is one of the easiest ways to render novel views whilst handling occlusions via z-buffering. The pipeline is this time straightforward (source to destination), since the world geometry (triangles) is directly projected and interpolated at the screen pixels (destination) yet the colours are interpolated in the reference texture. The main difference is that this requires a scattered interpolation which is computationally more expensive. The intensity warping can be defined such that

I * (P * ) = I Γ P * , E, w(V, T -1 , K * ) (8)
where E ∈ N J×3 contains J × 3 indices of each triangle (computed in the image as described in Section V-A.1), and Γ is the rasterisation function that interpolates the scattered vertices projected at the destination pixel locations P * . Rasterisation is directly implemented in hardware on GPUs and allows to efficiently perform scattered interpolation on a very large amount of triangles. With abuse of notation, the function

S * = S Γ P * , E, w(V, T -1 , K * ) (9)
will be a synthesis of the surface S at position T, where S * is the new predicted surface which contains new intensities, vertices, normals and uncertainties.

D. Warping comparisons

Forward warping (rasterisation) benefits from many hardware accelerated features that cannot be done efficiently with inverse warping (texture projection). In particular, mipmapping, anisotropic filtering (which prevents aliasing) and occlusions (via z-buffering). Even if the scattered interpolation can be highly optimised (see Section. V-A.1), the computational cost is less efficient than inverse warping, mainly due to triangulation and scattered interpolation of the surface.

On the other hand, inverse warping is only valid for small viewpoint changes where occlusions and perceptual aliasing can be neglected. In order to take advantage of hardware rasterisation, it is proposed here to use forward warping for view prediction and model update steps, which only have to be performed once per frame. Dense iterative registration, which requires multiple warps per frame, can use a fast inverse warping between the predicted frame and the live frame which are assumed to be close.

E. Multi-key-frame fusion

Supposing that an initial estimation of the current camera pose T is available (i.e. the last estimated pose), the surface prediction is performed by first extracting close key-frames from the graph. A simple criteria is used to select the M closest key-frames to the current frame based on the distance along the graph. This avoids choosing key-frames which are close geometrically but have not been connected visually during the acquisition (i.e. on the other side of a wall). Further connections could be made by performing loop closure and bundle adjustment on the entire key-frame graph.

Each key-frame is then rasterised and blended into a virtual frame at the predicted camera viewpoint such that

S * = M i=1 f S Γ P * , E, w(V i , T -1 T i , K * ) , (10)
where f (S(•)) is a blending function, that correctly fuses the synthesised images. In order to detect occlusions between the rendered surfaces before blending a Mahalanobis distance is computed between the candidate vertices to test their mutual dependency. Let us consider the sets {v * a , c * a } and {v * b , c * b } to be two vertices and uncertainties candidates for a new vertex v * . The distance is defined by

d M = e ⊤ 3 (v * a -v * b ) ⊤ (c * a + c * b ) -1 e ⊤ 3 (v * a -v * b ). ( 11
)
The blended vertex v * is obtained using a Chi-square test such that:

v * =      (v * a c * a -1 +v * b c * b -1 ) (c * a -1 +c * b -1 ) if d M < 3.841 v * a elseif e ⊤ 3 v * a < e ⊤ 3 v * b v * b else (12)
where 3.841 corresponds to 5% of error tolerance.

The intensities of the virtual image I * are also obtained by blending the reference key-frames. The weighting function proposed in [START_REF] Meilland | Super-resolution 3d tracking and mapping[END_REF] is used if the dependence test of [START_REF] Connolly | Cumulative generation of octree models from range data[END_REF] succeeds, otherwise the final intensity is assigned to the closest vertex. This function weights the relative resolution between the predicted image and the closest key-frames giving larger weights to images with closer viewpoints.

IV. MODEL UPDATE

This step consists in updating the model using the sensor pose estimation and the current image data. For volumetric approaches such as [START_REF] Newcombe | Kinectfusion: Real-time dense surface mapping and tracking[END_REF], this is achieved by integrating the new measurements (depths and intensities) into the TSDF through a weighted running average. For moving TSDF variants [START_REF] Whelan | Robust real-time visual odometry for dense RGB-D mapping[END_REF], [START_REF] Roth | Moving volume kinectfusion[END_REF], the volume can also be shifted if the camera pose exceeds some movement threshold.

For key-frame approaches such as [START_REF] Stückler | Model learning and real-time tracking using multi-resolution surfel maps[END_REF], new measurements can be integrated into nearby frames by warping the current data using the estimated pose, and a new key-frame is generated if the camera pose exceeds some movement threshold.

In the proposed approach, the current surface is rasterised onto the M key-frames that were used for the pose estimation. The integration proposed in [START_REF] Meilland | Super-resolution 3d tracking and mapping[END_REF] is employed to correctly fuse the depth-maps and intensities, allowing to reduce the noise contained in the raw depth images. A simple criterion for managing new key-frame generation is used in addition to a classic threshold on the camera motion. This criterion consists in monitoring the amount of occluded pixels between the current frame and the next frame prediction (which is computed according to the current pose estimate) using the test of equation [START_REF] Connolly | Cumulative generation of octree models from range data[END_REF]. A new key-frame is therefore added and connected to the M closest frames in the graph when the amount of occluded pixels in the image is greater than a defined threshold. This allows to dynamically scale the number of key-frames with the scene complexity.

V. EXPERIMENTS

A. Real-time implementation

The proposed approach was implemented and optimised for modern GPUs using the OpenCL library with OpenGL interoperability. This allows to use the OpenGL hardware pipeline for surface rasterisation with GLSL shaders and the more flexible OpenCL kernels for iterative registration with inverse warping. Extensive use of texture memory has been made to ensure fast random memory access with hardware bilinear interpolation. The iterative registration technique of Section II-B is performed with a coarse to fine multiresolution approach as detailed in [START_REF] Comport | Accurate quadrifocal tracking for robust 3d visual odometry[END_REF]. Since the RGB-D sensor usually provides noisy depth measurements, a bilateral filter is applied to remove noise whilst preserving discontinuities. The filtered depth-map is only used for pose estimation, while the raw depth-map is used for depth integration to preserve details in the integration process. The test platform used for experiments is a standard desktop PC running Ubuntu 12.10 with an Intel Core i7 2700k and a 2GB nVidia GeForce GTX670. Table I shows the average computing time of each step of the algorithm along with its standard deviation. The first step (initialisation) consists in loading the current RGB-D frame in GPU memory, performing bilateral filtering and computing the Gaussian pyramids of the images. The second step is the iterative pose estimation described in Section II-B. The third step is the model update of Section IV which consists in updating the key-frames. The last step is the next view prediction performed using the warping of Section III-C and the fusion of Section III-E. 1) Efficient triangulation on the GPU: Since the vertices of a key-frame are stored onto a pixel grid, one of the fastest way to rasterise a new surface is to use triangle strips. Triangle strips are a kind of primitive that specifies a series of connected triangles, sharing vertices, allowing for more efficient memory usage. However, such a surface does not take into account depth discontinuities, as well as missing values in the reference depth-map, which will create artefacts in the synthesis. One approach to handle this is to perform a simple 2D triangulation on the CPU using valid depth measurements and send the obtained indices to the GPU before drawing. This has some limitations:

• Since the depth-map is directly integrated in GPU memory, it must be sent to the CPU before triangulating; • The amount of triangle indices that have to be uploaded to the GPU after triangulation is generally enormous (up to 1843200 indices for a 640 × 480 To avoid performing this slow GPU/CPU memory transfer each time a reference view needs to be triangulated (i.e. multiple times at frame-rate), it is proposed here to perform the "triangulation" on the fly using geometry shaders. Geometry shaders are an optional part of the graphics rendering pipeline which is applied after primitive assembly. They allow to simultaneously access all the vertices of a primitive and can be used to discard or generate new primitives (only a few).

In this case, the entire set of reference vertices are drawn as points and sent to the vertex shader stage, which just sends the vertices positions to the geometry shader. The input of the geometry shader is then a vertex and the output is a triangle strip with a maximum of 4 vertices (2 triangles). The 3 neighbouring values of the current vertex are read via texture fetch and the triangle strip is emitted only if the following constraints are respected:

• All 4 depth values must be in a valid range (i.e. > 0).

• All triangle strip edge lengths must be smaller than a threshold to avoid depth discontinuities on the surface. The standard rasterisation pipeline is then applied with texture interpolation in screen space. Note that the rendered vertices and uncertainties are simultaneously interpolated and stored into RGBA floating point textures, using multiple render targets (MRT). Table II compares the computation time required for triangulating and rendering a 640 × 480 surface, containing 606404 triangles, using a naive CPU triangulation and the proposed approach. It can be seen that the full GPU triangulation is almost 14x faster than the naive CPU approach. Another advantage of triangulating on the fly is the memory footprint, since no temporary buffer of memory is used to store the triangles elements. Fig. 2. Bird eye view of a reconstructed office rendered in real-time using the proposed multi-key-frame fusion approach.

2) Reduction: To solve the linear system from equation [START_REF] Furukawa | Accurate, dense, and robust multiview stereopsis[END_REF], it is necessary to compute H = J T J and g = J T e. This is a large reduction problem since J is of dimensions 2nm × 6 and e is of dimensions 2nm × 1. Hence, H is of dimensions 6 × 6 and g is of dimensions 6 × 1. First note that the matrix H is symmetric, thus only 21 + 6 = 27 sums have to be computed for H and g. The final reduction can be interpreted as summing the 2nm rows of a 2nm × 27 matrix into a 1 × 27 vector that can be read back to the CPU before computing the pose update of equation [START_REF] Furukawa | Accurate, dense, and robust multiview stereopsis[END_REF]. Practically, in order save GPU cycles, only a partial reduction is performed on GPU using local memory, until a size of k×27, with k << 2nm. The final sum is then performed on the CPU and the system is inverted using Cholesky factorisation.

To further improve speed, SIMD instructions are used to perform 4 single precision instructions per cycle, reducing the problem from a 2nm × 27 to a 2nm × 7 reduction.

In the optimised implementation, the 2mn × 27 matrix is not explicitly constructed and stored in global memory. Instead the summation is performed in parallel as new error and Jacobian values are computed for each pixel. This allows to perform one iteration of the registration in only one kernel call without read/write operations via intermediate buffers.

B. Results

1) Quantitative results:

The proposed approach has been evaluated on the RGB-D benchmark of [START_REF] Sturm | A benchmark for the evaluation of rgb-d slam systems[END_REF]. Four selected sequences are reported on tables III, IV and V, which respectively represent the Absolute Trajectory Error (ATE), the Relative frame-to-frame Error (RPE) and the relative Root Mean Square Error per second (RMSE) corresponding to the drift per second. Three approaches have been compared: the Kintinuous volumetric approach of [START_REF] Whelan | Robust real-time visual odometry for dense RGB-D mapping[END_REF] with the ICP+RGB-D error metric, an image-based approach using only the Closest Key-frame in the Graph (CKF), and the proposed approach using a maximum of 5 View Predictions (VP). For the ATE and RPE, the median translation errors and the maximum translation errors are reported. It can be seen that for each error metric, both image-based approaches (VP and CKF) outperform the volumetric approach in all sequences. The proposed approach using frame prediction also performs better than the closest key-frame approach: reducing the amount of occlusions between the reference and the current view allows a better minimisation of the error in sensor space and reduces the maximum pose error and drift.

2) Qualitative results: The dense tracking and mapping algorithm has been successfully tested on a number of real scenes in real-time, ranging from small workspaces to large scale buildings. Figure 3 TOP shows a 3D reconstruction representing an entire floor of a building, reconstructed in real-time from a 100 meters trajectory using a hand-held RGB-D camera (Asus Xtion). The trajectory begins in the left alcove and the camera is moved across corridors and rooms, until the last alcove on the right. The environment contains many texture-less and geometrically unconstrained areas (corridors). Thanks to the bi-objective error minimisation, which combines a photometric minimisation and a geometric minimisation, the proposed approach is always able to robustly estimate an accurate pose. Figure 3 BOTTOM shows a side view of the final point cloud of Figure 3 TOP and demonstrates the very low drift of the approach since the entire reconstruction remains flat (floor). The final map, is only represented by 67 key-frames resulting in a point cloud of 20, 369, 566 3D points.

Figure 2 shows a synthesised global view of an office rendered using the proposed multi-key-frame 3D fusion approach using all the images of the graph. The rendered image is near photo-realistic which demonstrates the accuracy of the reconstruction pipeline. More detailed results can be found in the accompanying video (also available at: http://youtu.be/9p5IAkgh U4). In conclusion this paper has investigated two different representations commonly used for dense real-time localisation and mapping. The advantages and disadvantages of each approach have been demonstrated and a model has been proposed to unify these approaches. The proposed model is based on an efficient hybrid intensity and point-to-plane ICP formulation. Multiple reference images stored in a graph have been used to provide efficient occlusion handling to a multi-key-frame approach. In that respect a fast surface triangulation/rasterisation method has been developed on the GPU and it has been shown to use relatively low memory and be able map large scale buildings densely in real-time.

Future works will aim at detecting large loop closures using appearance-based techniques and correcting drift using pose-graph optimisation.

Fig. 3 .

 3 Fig. 3. Dense reconstruction of an entire floor obtained in real-time from a 100 meters trajectory containing 67 key-frames. The the graph of key-frames is superimposed on the images. TOP: a bird eye view of the reconstruction. BOTTOM: Side view of the reconstruction

TABLE I

 I

	ALGORITHM COMPUTING TIME
		Average time (ms)	Standard dev. (ms)
	Initialisation	12.90	2.49
	Pose estimation	11.78	2.06
	Map update/fusion	2.83	1.55
	View Prediction	3.14	1.38
	Total	30.63	3.08

TABLE II COMPARISON

 II OF THE COMPUTATIONAL TIME BETWEEN A NAIVE CPU

	TRIANGULATION AND THE PROPOSED GPU TRIANGULATION
		CPU+GPU (ms) Proposed (ms)
	Depth-map read-back	0.80	none
	2D Triangulation	10.50	none
	Elements update	3.56	none
	Rasterisation	0.77	1.13
	Total	15.64	1.13

TABLE III RELATIVE

 III ROOT MEAN SQUARE ERROR OF DRIFT IN METERS/SECOND (RMSE).

	Sequence	VP (m/s) CKF (m/s) Kintu. (m/s)
	fr1/desk	0.0259	0.0721	0.0393
	fr2/desk	0.0147	0.0161	0.0208
	fr1/room	0.0351	0.0502	0.0622
	fr2/large no loop	0.0695	0.1529	0.1795

TABLE IV ABSOLUTE

 IV TRAJECTORY ERROR IN METERS (ATE).

	Sequence	VP (proposed) Median Max Median Max Median Max CKF Kintu.
	fr1/desk	0.018 0.066 0.044 0.131 0.069 0.234
	fr2/desk	0.093 0.116 0.099 0.130 0.119 0.362
	fr1/room	0.144 0.339 0.201 0.430 0.158 0.421
	fr2/large no loop 0.187 0.317 0.228 0.437 0.256 0.878

TABLE V RELATIVE

 V FRAME-TO-FRAME POSE ERROR IN METERS (RPE). .0275 0.0046 0.1693 0.0045 0.0892 fr2/large no loop 0.0075 0.1620 0.0093 0.3324 0.0087 0.1101

	Sequence	VP (proposed) Median Max Median Max Median Max CKF Kintu.
	fr1/desk	0.0055 0.0390 0.0059 0.1899 0.0056 0.0655
	fr2/desk	0.0015 0.0098 0.0018 0.0171 0.0025 0.0108
	fr1/room	0.0041 0


Tommi Tykkälä
email: tommi.tykkala@lut.fi

Hannu Hartikainen
email: hannu.hartikainen@aalto.fi

Andrew I Comport
email: comport@i3s.unice.fr

Joni-Kristian Kämäräinen

RGB-D Tracking and Reconstruction for TV Broadcasts

Keywords: Dense Tracking, Dense 3D Reconstruction, Real-Time Tracking, RGB-D, Kinect

In this work, a real-time image-based camera tracking solution is developed for television broadcasting studio environments. An affordable vision-based system is proposed which can compete with expensive matchmoving systems. The system requires merely commodity hardware: a low cost RGB-D sensor and a standard laptop. The main contribution is avoiding time-evolving drift by tracking relative to a pre-recorded keyframe model. Camera tracking is defined as a registration problem between the current RGB-D measurement and the nearest keyframe. The keyframe poses contain only a small error and therefore the proposed method is virtually driftless. Camera tracking precision is compared to KinectFusion, which is a recent method for simultaneous camera tracking and 3D reconstruction. The proposed method is tested in a television broadcasting studio, where it demonstrates driftless and precise camera tracking in real-time.

INTRODUCTION

Rendering virtual elements, props and characters to live television broadcast combines augmented reality (AR) and video production. A robust camera estimation method is required for rendering the graphics from camera viewing direction. In film industry the process is known as matchmoving and it is traditionally done in post-processing [START_REF] Dobbert | Matchmoving: The Invisible Art of Camera Tracking[END_REF]. Matchmoving typically requires manual effort because the available tools are not fully automatic. In online broadcasting, special hardware based solutions exist for automatic camera tracking, but their value is degraded by limited operating volume, weaker reality experience and expensive price.

The goal of this work is to develop an affordable, portable and easy-to-use solution for television production studios. Marker-based AR techniques, such as ARToolKit [START_REF] Kato | Marker tracking and hmd calibration for a video-based augmented reality conferencing system[END_REF], are trivial to use but visible markers are irritating in the studio scene. Thus, we seek solution from more recent techniques which are able to track the camera without any markers. Visual simultaneous localisation and mapping (vSLAM) techniques are a true option, because they track the camera using a 3D model which is concurrently built based on visual measurements [START_REF] Davison | MonoSLAM: Real-time single camera SLAM[END_REF]. Recently, low-cost RGB-D sensors have developed to a level where they can provide realtime stream of dense RGB-D measurements which are directly usable for camera tracking and scene reconstruction purposes.

In this work, we introduce an image registration based dense tracking and reconstruction method particularly for TV studios. Tracking precision is good, but over long sequences, time-evolving drift will eventually displace virtual props. We avoid drift by generating a RGB-D keyframe model, and tracking the camera relative to the nearest keyframe. The results are demonstrated with real studio shots. We thank Heikki Ortamo and Jori Pölkki for their professional support in a TV production studio.

Related Work

Traditionally vSLAM methods detect and track a sparse set of feature points which are matched in several images. The first feature-based visual SLAM methods used the extended Kalman filter (EKF) to update pose and structure [START_REF] Davison | MonoSLAM: Real-time single camera SLAM[END_REF]. However, bundle adjustment has replaced EKF, because it is more accurate. PTAM [START_REF] Klein | Parallel tracking and mapping for small ar workspaces[END_REF] separated tracking and mapping into two parallel modules, where the mapping part was essentially bundle adjustment. To avoid feature extraction and matching problems, the raw pixel measurements can be used directly. DTAM was introduced as the dense version of PTAM, which allows every pixel to contribute to pose and structure estimation (Newcombe et al., 2011a). Finally, with KinectFusion system which replaces monocular camera by a RGB-D sensor [START_REF] Newcombe | Kinectfusion: Realtime dense surface mapping and tracking[END_REF], pose tracking and structure estimation have become mature enough to be considered for live TV broadcasts.

In this work, we adopt the dense RGB-D approach [START_REF] Tykkala | Direct iterative closest point for real-time visual odometry[END_REF][START_REF] Audras | Real-time dense rgb-d localisation and mapping[END_REF][START_REF] Comport | Accurate quadri-focal tracking for robust 3d visual odometry[END_REF]. Our method differs rather strongly from the KinectFusion, because we estimate camera pose using dense RGB-D measurements instead of depth maps only. This is necessary because studio settings often contain planar surfaces for which the KinectFusion fails to track the camera. We do not estimate 3D structure concurrently because it can be solved prior to broadcasting. This simplification lightens computation and enables running the system with a low-end hardware. Robustness to outlier points is increased by selecting the photometrically stable pixels and by omitting unreliable regions (occlusion and moving objects) with the M-estimator.

DENSE TRACKING METHOD

Cost Function

Our dense pose estimation is defined as a direct color image registration task between the current image I : R 2 ⇒ R and a frame K * = {P * , c * }, where P * = {P 1 , P 2 , . . . , P n } is a set of 3D points and c * = {c 1 , c 2 , . . . , c n } are the corresponding color intensities. Our goal is to find the correct pose increment T(ω, υ) which minimizes the following residual:

e = I (w (P * ; T(ω, υ))) -c * , (1)
where w(P ; T) is a projective warping function which transforms and projects P into a new view using the 4 × 4 transformation matrix T and the intrinsic matrix K (constant, omitted in the notation). The exact formula is

w(P ; T, K) = w P ; R t 0 1 , K = N(K(RP + t)),
(2) where N(p) = (p 1 /p 3 , p 2 /p 3 ) dehomogenizes a point. T(ω, υ) is defined as the exponential mapping which forms a Lie group T(ω, υ) ∈ SE(3) [START_REF] Ma | An invitation to 3-D vision: from images to geometric models[END_REF]:

T(ω, υ) = e A(ω,υ) , A(ω, υ) = [ω] × υ 0 0 , (3)
where ω and υ are 3-vectors defining rotation and translation. From the practical point of view it is convenient to re-define (3) as

T(ω, υ) = Te A(ω,υ) , (4)
which generates smooth increments around the base transform T. This allows using iterative optimization and particularly the inverse compositional trick for estimating the transformation efficiently [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF].

Minimization

We adopt the inverse compositional approach for efficient minimization of the cost function [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF]. The cost is reformulated as

c * (ω, υ) = I * w P * ; e -A(ω,υ) e = c * (ω, υ) -I w P * ; T ,
where reference colors c * are now a function of the inverse motion increment. The benefit is gained when computing the Jacobian

J i j = ∆I * (w(P i ; I)) ∂w(P i ; I) ∂x j , (5)
where

∆I (p) = [ ∂I (p) ∂x ∂I (p) ∂y ]
. By this trick J does not depend on the current transformation anymore and it is sufficient to computed it only once for each K * .

Gauss-Newton iteration is well-suited for minimization when the initial guess is near the minimum and the cost function is smooth. The cost function is locally approximated by the first-order Taylor expansion e(x) = e(0) + Jx, where e(0) is the current residual. Now the scalar error function becomes

1 2 e(x) T e(x) = 1 2 e(0) T e(0) + x T J T e(0) + 1 2 x T J T Jx (6)
where the derivative is zero when

J T Jx = -J T e(0) . (7)
Because J T J is a 6 × 6 positive definite matrix, the inversion can be done efficiently using Cholesky method. The matrix multiplication associativity in SE(3) enables collecting the previous increments into the base transform by

T k • • • T 1 e A(ω,υ) ⇒ Te A(ω,υ) (4).

Keyframe Based Reference K *

Minimisation of the cost function (1) with the Gauss-Newton method in Sec 2.2 provides us real-time dense RGB-D camera tracking (and scene reconstruction), which frame by frame finds the optimal transformation T(ω, υ) between a previous frame K * and a current frame (image I ). This approach has one disadvantage: small per frame estimation error cumulates into global drift which gradually displaces virtual elements. Fortunately, the studio scene can be prerecorded into a set of keyframes prior to broadcasting (see Fig. 1). The tracking is then defined relative to the nearest keyframe K * ← select closest(K * j ). For small studios, a single dense tracking sweep is sufficient, but methods exist also for building larger models [START_REF] Henry | RGB-D mapping: Using Kinect-style depth cameras for dense 3D modeling of indoor environments[END_REF].

Robust estimation

Now we know that our reference keyframes contain merely static geometry with Lambertian reflection, it is possible to increase robustness of the estimation by emphasising stronger color and depth correlation between the static reference points and the current 3D points matched in each iteration. We do this by employing M-estimation in which uncertainty based weights w k ∈ [0, 1] ∈ R are given to the residual elements e k . Color correlation is enforced by the Tukey weighting function

u k = |e k | c * median(e a ) , (8)
w c k = (1 -( u k b ) 2 ) 2 if |u k | <= b 0 if |u k | > b , (9)
where e a = { e 1 , . . . , e n }, c = 1.4826 for robust standard deviation, and b = 4.6851 is the Tukey specific constant.

The depth correlation weights are computed simultaneously by a depth lookup

e z = Z (w(P ; T(ω, υ))) -e T 3 T(ω, υ) P 1 , (10)
where Z : R 2 ⇒ R is the depth map function of the current RGB view and e T 3 = (0, 0, 1, 0) selects the depth coordinate. When the standard deviation of depth measurements is τ, the warped points whose depth differs more than τ from the current depth map value can be interpreted as foreground actors/outliers. Thus we define a diagonal matrix

w z k = max (1 -e 2 z (k)/τ 2 ), 0 2 , ( 11
)
The weighted Gauss-Newton step is obtained by re-writing (7)

J T WJx = -J T We, ( 12
)
where W is a diagonal matrix with diag(W) k = w k c w k z . The robust step is obtained by J ⇐ √ WJ and e ⇐ √ We. The both weight components are quadratic and therefore the square roots are not evaluated in practise.

Selecting keyframe points

The reference points P * can be freely selected from the keyframe. When considering the minimization of the photometrical error, only the 3D points with image gradients infer pose parameters. Therefore the majority of the points can be neglected based on the magnitude of the image gradient. Gradient vector ∇c is evaluated at the projection of P * using bilinear interpolation. Instead of sorting |∇c|, we generate the histogram of the magnitudes. The value range is bounded to [0, 255]. We seek the bin B t for which Σ 255 k=B t -1 h(k) < n <= Σ 255 k=B t h(k), where n is the number of points to be selected.

IMPLEMENTATION DETAILS

Our method was implemented in Ubuntu Linux environment using open software tools, Kinect RGB-D sensor and a commodity PC laptop hardware on which the method runs real-time. The main limitations are the operating range (≈ 1m -5m) and that it uses controlled IR lighting which may not work in outdoors. Below, we discuss Kinect related implementation issues.

Kinect Calibration

Since Microsoft Kinect can be set into a special mode where the raw IR images can be stored, it is possible to use standard stereo calibration procedure for obtaining the calibration parameters for the IR and the RGB view [START_REF] Bouguet | Camera calibration toolbox for matlab[END_REF] (Fig. 2). IR view and the depth view are trivially associated with image offset (-4, -3).

A single camera calibration is used to initialize IR and RGB camera parameters. Then calibration is followed by a stereo procedure which re-estimates all free parameters (K IR , K RGB , kc RGB ,T b ). The lens distortion parameters of the IR camera are forced to zero, because data has already been used to generate the raw disparity map. This means that the IR lens distortion is compensated by tweaking other parameters. The RGB camera lens distortion parameters kc RGB are estimated without any special concerns. In practice, the distortion seems to be minor and the first two radial coefficients are sufficient (kc RGB = (0.2370, -0.4508, 0, 0, 0), and the stereo baseline is b = 25.005mm. T b stores the baseline transform as 4×4 matrix. The conversion from raw disparities into depth values can be done by z = 8p f B-d , where p is the baseline between the projector and the IR camera, B is a device specific constant and f is IR camera focal length in pixel units. Note, that this reconstruction method is merely an approximation which procludes measurements at long ranges. There are also dedicated calibration toolboxes for RGB-D sensors, which model disparity distortion accurately [START_REF] Herrera | Joint depth and color camera calibration with distortion correction[END_REF].

Dense Tracking with Kinect

The The cost function is minimized using coarse-tofine approach using image pyramid with 80 × 60, 160×120 and 320×240 layers for each RGB-D input frame.

EXAMPLES

Dense tracking vs Kinfu

Kinfu is the open source implementation of Kinect-Fusion [START_REF] Rusu | 3D is here: Point Cloud Library (PCL)[END_REF]. We compare our dense tracking accuracy with Kinfu using the RGB-D SLAM benchmark provided by Technical University of Münich [START_REF] Sturm | Towards a benchmark for rgb-d slam evaluation[END_REF] (Figure 3). Dense tracking is executed incrementally by using a recent view as the reference. Thus, both methods aim at tracking the camera pose without a prior model and small drift will be present. The major difference between the systems is that Kinfu optimizes a voxel based 3D structure online and uses the iterative closest point (ICP) approach for pose estimation. Kinfu has small drift when the voxel size is small and geometry is versatile. Kinfu fails in bigger operating volumes, because voxel discretization becomes coarse and, for these sequences especially, the volume will also contain planar floor which break downs ICP 1 . In broadcasting studios, scenes are often larger than (3m) 3 and geometrical variations can not be guaranteed. Our dense tracking demonstrates robust tracking even when planar surfaces are present, because the cost function matches also scene texturing. Memory consumption is low even in larger operating volumes, because RGB-D keyframes can be memoryoptimized based on the viewing zone.

Table 1 shows the comparison between our method and Kinfu numerically.

The dense tracking drifts 1.08cm/s with the slower freiburg2/desk sequence and 2.60cm/s with the faster freiburg1/desk sequence. Our dense tracking has smaller error when the camera is moving faster. Kinfu has lower drift when using (3m) 3 voxel grid, but fails to operate in bigger volumes (8m) 3 that match with broadcasting studios. Also the computational requirements of Kinfu are significantly higher compared to our approach even though Kinfu is executed on a powerful desktop 1 Video: http://youtu.be/tNz1p1sdTrE Table 1: Drift is evaluated using standard sequences with known ground truth trajectories. The proposed approach has smaller drift when a camera is moving faster. The keyframe model can be memory-efficiently built. Kinfu has smaller drift when a camera is moving slower but bigger scenes are not possible due to inscalability of the voxel grid. Problems exist with planar surfaces. The computational requirements are higher even though Kinfu is executed on a powerful desktop hardware.

hardware. Drift was measured by dividing the input frames into subsegments of several seconds (10 and 2 correspondingly) whose median error was measured against the ground truth. 1 second average error was computed from the median subsegment. The error values are computed from bigger windows to average out random perturbations and neglect gross tracking failures which occur with Kinfu in all cases except on freiburg2/desk using (3m)3 volume.

Driftless keyframe tracking

The relative poses between the keyframes could, in theory, be obtained by bundle adjustment techniques [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] if feature point extraction and matching succeeded, and a good initial guess would exist. With the studio sequences recorded, texturing was so limited that popular tools, such as Bundler [START_REF] Snavely | Photo tourism: Exploring photo collections in 3d[END_REF] failed without manual annotated feature points. As an alternative solution, we utilised the proposed dense RGB-D tracker to incrementally build a keyframe model. It is noteworthy, that a single quick sweep of the scene produces accurate keyframe model2 . Figure 4 illustrates how dense tracking is used to sweep a keyframe model of the studio scene rapidly. With pre-recorded keyframes, the online broadcasts are guaranteed to operate with the correct 3D geometry. K * j were selected by picking RGB-D frames evenly from the recorded sequence. A keyframe database is illustrated in Figure 1.

In Figure 5, we show how the drift increases in a studio environment when using dense tracking. On the right-hand side, the drift problem is solved by tracking relative to keyframes 3 . Note that in longterm use, even a small number of keyframes eventually outperforms the dense tracking due to drift. However, a small keyframe number produces drift jumps between the keyframes which can be visually disturb-

CONCLUSIONS

In this work, camera pose tracking was defined as a photometrical registration problem between a reference frame and the current frame. To remove the global drift in incremental tracking, the closest prerecorded keyframe was chosen to be the motion reference. The system was designed to be an affordable solution for TV broadcasting studios relying only on the Kinect sensor and a commodity laptop. The proposed approach performs robustly in a standard benchmark, where KinectFusion has problems with planar surfaces and limited voxel grid resolution. Our future work will address the practical issues how studio staff and camera men can use our computer vision system in live broadcasts. Moreover, combination of the best properties of our approach and KinectFusion will be investigated.

Figure 1 :

 1 Figure 1: Stored scene keyframes and their pixels mapped to a common 3D world model.

  p and B are estimated by solving the linear equation -1 Z A B T = D, where A = 8p f , Z is n × 1 matrix of reference depth values z k from the chessboard pattern (Caltech calibration), and D is a n × 1 matrix whose elements are d k z k . The parameters will be p ≈ 75mm and B ≈ 1090.

Figure 2 :

 2 Figure 2: RGB and IR images of the calibration pattern.

  reference point clouds {P * , c * } in (1) were generated from Kinect RGB image and disparity map in the following way. Kinect Bayer images were converted into RGB format, downsampled into 320x240 and undistorted from lens distortions. Downsampling is almost lossless due to sparsity of RGB values in the Bayer pattern. The raw disparity map was first converted into a depth map, downsampled into 320x240 size using max filter and then transformed into a point cloud P * IR . Maximum filtering is chosen because it does not produce artificial geometry. P * was then generated from T b P * IR , where T b is the baseline transformation between the IR and RGB cameras. Points p ∈ P * do not exactly project to the pixel centers of the RGB image grid, and thus, bi-linear interpolation is used for generating the corresponding intensities c * .

Figure 3 :

 3 Figure 3: On the left Kinfu result for freiburg desk2 sequence using 3x3x3, 5x5x5 and 8x8x8 meter voxel volume. Kinfu gains lower drift due to structural optimization, but planar surfaces cause tracking failures. 3x3x3 volume does not contain floor and therefore Kinfu works well. Limited operating volume is a problem for practical use in the studio. On the right, our dense tracking drift illustrated when a fixed keyframes are not used. Problems with planar surfaces or limited operating volume do not exist. Green dots reveal the selected points for given RGB-D measurement.

Figure 4 :

 4 Figure 4: Camera trajectory solved from Kinect input by minimizing the cost function (1). A small shark is rendered into the studio scene.

Figure 5 :

 5 Figure5: Camera moving front and back along a fixed 3.30m studio rail. On the left, three images taken from the beginning, middle and the end of the rail. Green regions illustrated selected points. On the right, comparison between dense tracking and keyframe tracking. In dense tracking drift increases in time, but keyframe tracking maintains small bounded error. A person is moving in the scene during the last four cycles.

Video: http://youtu.be/wALQB3eDbUg

Video: http://youtu.be/zfKdZSkG4LU


Patrick Maget
email: patrick.maget@cea.fr

Olivier Février

Xavier Garbet

Hinrich Lütjens

Jean-François Luciani

Alain Marx

Jean-Francois Luciani

Corrigendum : Extended Magneto-Hydro-Dynamic model for Neoclassical Tearing Mode computations

An extended fluid model covering neoclassical physics has been implemented in the XTOR code for Magneto-Hydro-Dynamic computations and is described in [START_REF] Maget | Extended magneto-hydro-dynamic model for neoclassical tearing mode computations[END_REF]. This model allows recovering neoclassical flux average quantities like the bootstrap current and the poloidal ion flow, and in the presence of a magnetic island, it generates a drive for Neoclassical Tearing Modes. The contribution of parallel heat fluxes on the bootstrap current is significant, and it was retained in the simulations presented in this paper. However, we have realized that the closure that is used for these parallel heat fluxes does not meet an important constraint on its spatial distribution, and although this does not change the equilibrium quantities, we do see an impact on the dynamics of the magnetic island. In the following, we propose a different closure that satisfies this constraint, and we also present a slight modification of the neoclassical implementation in the momentum equation. In order to show how these modifications impact the simulations, we present a limited number of examples performed with the new model, with diagnostics tools allowing a better understanding of the physics at play.

Corrections to the MHD model

Closure on parallel heat fluxes

The modifications that we report here concerns the expression of parallel heat flows and the momentum equation. The parallel heat flux of the species a can be expressed as a function of plasma flows and neoclassical forces from each species b. The flux average parallel heat flux satisfies :

u 2 a = b C ab u 1 b + C ab 1 S b 1 + C ab 2 S b 2 (1)
with

u 1 b ≡ V b • B and u 2 a ≡ 2q a • B/(5p a )
. The expressions for C ab , C ab 1 and C ab 2 are given in the appendix of reference [START_REF] Maget | Extended magneto-hydro-dynamic model for neoclassical tearing mode computations[END_REF]. The terms S b 1,2 need however some modification. Indeed, the expressions given in reference [START_REF] Maget | Extended magneto-hydro-dynamic model for neoclassical tearing mode computations[END_REF] have a dependence in sin 2 θ with θ the poloidal angle, due to the poloidal variation of the magnetic field amplitude as cos θ.

Since we have no evolution equation for the heat flux u 2 a , this strong poloidal variation cannot be regularized and leads to unphysical poloidal distribution of the parallel heat flux. In order to solve this difficulty while maintaining the correct flux averaged relation that leads to the equilibrium bootstrap current, we adopt the following definitions :

S a 1 S a 2 = B 2 B • ∇θ V * a + V E V * Ta • ∇θ (2)
with V * T a = B × ∇T a /(e a n a B 2 ), with a perpendicular heat flux given by 2q a⊥ /(5p a ) = V *

Ta

(3)

The closure on the parallel heat flux is then derived from equilibrium constraints. It can be written as

u 2 a = 2 5p a q a • ∇θ B • ∇θ B 2 -S a 2 (4)
At equilibrium, assuming that the heat flux is divergence free and has no radial component, we have

∂ θ q a • ∇θ B • ∇θ = 0 (5)
so that this quantity is a flux function. It is also easily shown that S a 2 is a flux function at equilibrium, so that we can write

u 2 a = u 2 a B 2 B 2 + B 2 B 2 -1 S a 2 (6)
We will use this relation as the closure on the parallel heat flux. The contribution of the heat flux in the pressure anisotropy (see equation (29) of [START_REF] Maget | Extended magneto-hydro-dynamic model for neoclassical tearing mode computations[END_REF]) is therefore :

V * T a + u 2 a B 2 B • ∇θ = B 2 B 2 V * T a • ∇θ + u 2 a B 2 B • ∇θ (7)
This formulation gives therefore the correct contribution of the flux averaged parallel heat fluxes but ignores how these fluxes can participate in the neoclassical island drive.

Momentum equation

Regarding the ion neoclassical friction force in the momentum equation, we used the following expression in [START_REF] Maget | Extended magneto-hydro-dynamic model for neoclassical tearing mode computations[END_REF] :

∇p + ∇ • (Π i + Π e ) = p -p ⊥ [(∇ • b) b + κ] + ∇ p -p ⊥ + ∇p ⊥ (8)
with p = p i + p e . This formulation remains satisfactory for moderate anisotropy where p ⊥ can be identified with p. However, this modification has proved to be unnecessary from a computation performance point of view, so that we now adopt the expression

∇•(Π i + Π e ) = p -p ⊥ [(∇ • b) b + κ]+∇ p -p ⊥ - 1 3 ∇ p -p ⊥ (9)
that allows better comparison with formal derivations of neoclassical island theory (as in [START_REF] Smolyakov | On neoclassical effects in the theory of magnetic islands[END_REF] for example).

Measure of neoclassical contribution

In order to characterize the neoclassical drive in a simple way, we define a scalar R(f, g) that quantify the correlation between two 3D fields f and g over the plasma volume V :

R(f, g) = dV f g / dV f 2 1/2 dV g 2 1/2 (10)
This scalar is unity if f is proportional to g with the same sign, and minus unity if the sign is opposite. If the two fields are completely uncorrelated, R(f, g) should be close to zero. The link between the bootstrap current perturbation and the pressure gradient perturbation can then be evaluated by taking f ≡ [J bs,ϕ ] n =0 and g(h) ≡ [-∂ y h] n =0 , using h = p (pressure), T (temperature), or ρ (density), y ≡ √ ψ the radial co-ordinate and n the toroidal mode number. In the following, we note R bs (h) = R(f, g(h)). In order to quantify the contribution of the perturbed E × B velocity V E , we also compute

R bs (v) = R(f, g) with g = [V E,θ ] n =0 .
Although the bootstrap current can be formally expressed as a sum of the contributions from g(p), g(ρ) and g(v), note that the sum of their correlation terms R bs is not constrained to be unity.

In the evaluation of the bootstrap drive to island growth, it is generally assumed that the neoclassical friction is strong enough to force the perturbed flow to its neoclassical drive. This motivates the use of an ad-hoc bootstrap model where J bs ∝ -∇p for non linear simulations of NTM dynamics [START_REF] Lütjens | Nonlinear three-dimensional MHD simulations of tearing modes in tokamak plasmas[END_REF][START_REF] Lütjens | Linear and nonlinear thresholds of neoclassical tearing modes in tokamaks[END_REF][START_REF] Popov | Simulation of neoclassical tearing modes (NTMs) in the DIII-D tokamak. I. NTM excitation[END_REF][START_REF] Maget | Modelling of (2,1) NTM threshold in JET advanced scenarios[END_REF]. With such a model, we have naturally R bs (p) ≡ 1 and a bootstrap perturbation that is closely linked to the island structure, due to the large parallel diffusivity that equilibrates the temperature along magnetic field lines.

Numerical simulations

We have performed numerical simulations on the same equilibrium that was used in [START_REF] Maget | Extended magneto-hydro-dynamic model for neoclassical tearing mode computations[END_REF], with an ion density n i (0) = 2 × 10 19 m -3 , electron temperature T e (0) = 3910eV , and τ ≡ T i /T e = 1, a Lundquist number S 0 = 10 7 , a magnetic Prandtl number Prm≡ ν/η = 10, χ /χ ⊥ = 10 8 , χ ⊥ /η = 150, and D ⊥ = 2/3χ ⊥ .

In oder to save computation time, a (2,1) seed island is inserted at the beginning of the simulation, with a size W = 4%. The evolution of the island width, the bootstrap current fraction at the resonance and the correlations R bs (p) and R bs (v) are shown in figure 1 for different amplifications of the bootstrap current. The simulation for f bs = 1 can be compared with figure 7 of [START_REF] Maget | Extended magneto-hydro-dynamic model for neoclassical tearing mode computations[END_REF]. Although the saturation size is similar, we note the absence of the island breaking events, similar to the sporadic growth of plasmoids, that were present with the previous closure. The scan in bootstrap current shows, as noted before, a weak effect that is even opposite to the usual destabilizing drive. This can be understood with the correlation between the bootstrap perturbation and the island structure R bs (p) that is slightly negative. On the contrary, the correlation R bs (v) is slightly positive, showing that the contribution of the perturbed electrostatic potential in the bootstrap current is larger than that of the diamagnetic flows.

This result can be understood by the fact that temperature flattening inside the saturated island is small, thus limiting the nonlinear impact of the island on the bootstrap current. In the simulation, this temperature flattening is controlled by the characteristic transport width W χ [START_REF] Fitzpatrick | Helical temperature perturbations associated with tearing modes in tokamak plasmas[END_REF] defined in [START_REF] Maget | Extended magneto-hydro-dynamic model for neoclassical tearing mode computations[END_REF]. We have verified that the correlation R bs (p) increases when the characteristic width W χ is reduced, as shown in figure 2 where it becomes slightly positive. Whereas W χ was originally larger than the saturated island width (W χ = 4.7%), it has been decreased to W χ = 3.3% by increasing the parallel diffusivity coefficient by a factor of 4. At the same time, the curvature stabilization effect is reinforced and reduces the island saturation [START_REF] Lütjens | Curvature effects on the dynamics of tearing modes in tokamaks[END_REF].

But the insufficient temperature flattening is not the only mechanism causing the vanishing of the bootstrap drive that is observed in the simulation of figure 1. A strong contributor is found to be diamagnetic rotation, that opposes diffusive processes in the equilibration of temperature along field lines. For changing the level of diamagnetic rotations without modifying the magnetic equilibrium, we increase plasma density, and decrease temperature accordingly to keep the same pressure. By doing so, we also modify the bootstrap current and the neoclassical friction coefficients (see table 1). The amplitude of this friction term is expected to be crucial for the neoclassical drive following [START_REF] S V Konovalov | Role of anomalous transport in onset and evolution of neoclassical tearing modes[END_REF][START_REF] Mikhailovskii | An analytic approach to developing transport threshold models of neoclassical tearing modes in tokamaks[END_REF] where a characteristic width W ν ∝ (ν/µ i ) 1/2 was introduced. With our notations and normalizations, this width is given by:

W real ν = 1 2 (ν/µ i ) 1/2 = 1 2 Prm 1/2 S real µ i -1/2 (11)
In our simulations, we do not rescale the ion neoclassical friction term µ i with the Lundquist number, but we have the possibility to scale it by a factor α µ . We have therefore

W sim ν = S real S 0 α µ 1/2 W real ν (12)
The values for W real ν and W sim ν corresponding to the simulations are reported in table 1. While in real situations the visco-neoclassical width is weakly sensitive to the scan in density that we have performed, a wide range is covered in the simulations, mainly due to the variation of the real Lundquist number. In particular, the ion neoclassical friction coefficient is increased by one order of magnitude between the cases N i (0) = 2×10 19 m -3 and N i (0) = 8 × 10 19 m -3 . By setting α µ = 10 (or even α µ = 20) in the reference case at N i (0) = 2 × 10 19 m -3 , we reduce the visco-neoclassical width below the island width, but we remain still larger than the real one. In this range at least, the correlation parameter remains constant and negative (see figure 4), showing that the dominant effect seen in the simulations comes from the change of diamagnetic rotations. At a lower level of diamagnetic rotation, the pressure contribution to the bootstrap drive increases significantly (see figure 3). This is confirmed by the fact that the island saturation increases when the bootstrap current is amplified by the factor f bs (see figure 1, right plot).

Case

S real N i (0) Lundquist number consistent with plasma parameters, central ion density, ratio of the equilibrium bootstrap current to the local current density at q = 2, visco-neoclassical width consistent with plasma parameters and the one used in the simulations. The series of simulations shown here demonstrate that we have now a neoclassical model that recovers qualitative expectations from theory, i.e. a bootstrap drive that increases the island saturation when it is sensitive to the island structure (positive R bs (p)). The importance of the electrostatic contribution to the bootstrap current has been clearly evidenced when diamagnetic rotations are large enough, leading to a strong reduction of the bootstrap drive. The role of the visco-neoclassical width W ν was not evidenced in the situation that we have encountered here, but it might be important in cases where the bootstrap drive is stronger.

J bs /J 0 W real ν W sim ν [×10 7 ] [
In addition to the study of neoclassical tearing modes, this model also allows in a consistent approach to investigate the transport of impurities, following the generic mechanism involving the parallel friction force between the main ions and trace species [START_REF] Ahn | Dynamics of heavy impurities in non-linear MHD simulations of sawtoothing tokamak plasmas[END_REF].

Figure 1 .

 1 Figure 1. Influence of the bootstrap current on the island dynamics at low density (left, reference case), and high density (right), with f bs the amplification factor. From top to bottom: island size, bootstrap current fraction at q = 2, correlation indicators R bs (p) and R bs (v).

Figure 2 .

 2 Figure 2. Influence of the characteristic transport width W χ on the island dynamics at low density (left, reference case), and high density (right), with W χ = 4.7% for χ /χ 0 ⊥ = 10 8 and W χ = 3.3% for χ /χ 0 ⊥ = 4 × 10 8 (χ 0⊥ is the value of χ ⊥ at plasma center, and it is about a factor two above at q = 2). From top to bottom: island size, correlation indicators R bs (p) and R bs (v).

Figure 3 .

 3 Figure 3. Influence of plasma density on the island dynamics. From top to bottom: island size, bootstrap current fraction at q = 2, correlation indicators R bs (p) and R bs (v).

Figure 4 .

 4 Figure 4. Influence of the visco-neoclassical width W ν on the island dynamics for the reference case. We have W ν = 14.0% for α µ = 1 and W ν = 4.4% for α µ = 10. From top to bottom: island size, correlation indicators R bs (p) and R bs (v).

Table 1 .

 1 10 19 m -3 ] at q = 2 Characteristic parameters (Prm=10, S 0 = 10 7 for the simulations):

	Ref	75.8	2	0.22	0.016 0.140
	Ref f bs = 0.5	75.8	2	0.11	0.016 0.140
	Ref f bs = 1.5	75.8	2	0.33	0.016 0.140
	Ref, α µ = 10	75.8	2	0.22	0.016 0.044
	Ref, α µ = 20	75.8	2	0.22	0.016 0.031
	N ref × 2	20.2	4	0.20	0.015 0.066
	N ref × 3	9.3	6	0.18	0.016 0.047
	N ref × 4	5.4	8	0.15	0.017 0.041
	N ref × 4 f bs = 0.5	5.4	8	0.075	0.017 0.041
	N ref × 4 f bs = 1.5	5.4	8	0.225	0.017 0.041

Acknowledgments

This work has been carried out within the framework of the French Research Federation for Fusion Studies. We benefited from HPC resources from GENCI (project 056348) and from and from Marconi Aix-Marseille Université project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements dAvenir" supervised by the Agence Nationale pour la Recherche. This work is part of the project AMICI funded by the Agence Nationale pour la Recherche (ANR-14-CE32-0004-01).