Tommi Tykkälä

Andrew I Comport

Joni-Kristian Kämäräinen

Hannu Hartikainen

Live RGB-D Camera Tracking for Television Production Studios

Keywords: Dense Tracking, Dense 3D Reconstruction, Augmented Reality, Real-Time, RGB-D, GPU

In this work, a real-time image-based camera tracker is designed for live television production studios. The major concern is to decrease camera tracking expenses by an affordable vision-based approach. First, a dense keyframe model of the static studio scene is generated using image-based dense tracking and bundle adjustment. Online camera tracking is then defined as registration problem between the current RGB-D measurement and the nearest keyframe. With accurate keyframe poses, our camera tracking becomes virtually driftless. The static model is also used to avoid moving actors in the scene. Processing dense RGB-D measurements requires special attention when aiming for real-time performance at 30Hz. We derive a real-time tracker from our cost function for a low-end GPU. The system requires merely a RGB-D sensor, laptop and a low-end GPU. Camera tracking properties are compared with KinectFusion. Our solution demonstrates robust and driftless real-time camera tracking in a television production studio environment.

Introduction

In this work, we develop a real-time image-based camera tracking solution for television broadcasting studios. Such a system enables adding virtual scene elements and interactive characters into live television broadcasts. The film industry knows this technique as matchmoving, which is traditionally done in postprocessing using semi-automatic trackers [START_REF] Dobbert | Matchmoving: The Invisible Art of Camera Tracking[END_REF]. A tracker is a tool which estimates a 3D camera motion trajectory based on observed 2D point tracks [START_REF] Laboratorium | für Informationstechnologie, Hannover, Voodoo camera tracker: A tool for the integration of virtual and real scenes[END_REF]. A tracker should only use the static scene geometry in order to avoid estimation bias. Scene segmentation is not typically fully automatic process and studio-specific standard technologies such as bluescreens, image masks and movement zones are used. Online matchmoving solutions exist, such as MotionAnalysis motion capture system, which tracks a camera in real-time based on passive markers attached to the camera. Brainstorm Multimedia product is compatible with many camera trackers and can be used to render view-dependent graphics in real-time. Although motion capture is precise and enables a large operating volume, the total price of such a system is currently 200 -500k€. Overall, matchmoving is a time-consuming process, which is expensive and lacks automatic, affordable, and easy-to-use tools.

The main aim in this paper is to reduce the costs of matchmoving in studio environments by introducing an affordable RGB-D sensor based camera tracking tool which operates in real-time, fits studio use, and merely requires a low-end GPU. Our contributions are the following. 1) A RGB-D keyframe-based tracking method is proposed, which does not suffer from time-evolving global drift. A static studio scene is first modeled as a database of RGB-D keyframes, which are obtained by using our dense tracking approach and bundle adjustment. The database is then used as a reference for realtime pose estimation (Figure 1) 1 . By defining the camera tracking problem relative to the nearest keyframe, the time-evolving drift present in various other systems is avoided. We show how the keyframe tracking eventually outperforms incremental dense tracking. The estimation is robust due to gradient-based pixel selection and a M-estimator.

When aiming at processing dense RGB-D data in real-time, the algorithms must be parallelized to obtain sufficient performance and scalability properties. 2) The Figure 1: Top: A RGB sensor is used to build a keyframe database which contains the views in the camera motion zone. The TV camera pose is estimated by registering the image content with the nearest keyframe. Bottom: A real-time AR broadcast with virtual 3D characters which are rendered from the estimated pose. See example video 1 .

computational scalability and performance of the camera tracking is improved by designing the full algorithm for a commodity GPU. A detailed description from our cost function definition into an efficient GPU implementation is given. 3) With a static keyframe-based 3D model available, the dynamic foreground points can be rejected from the camera pose estimation by observing discrepancies in intensity and depth simultaneously. Generally in RGB-D tracking, outliers both in color and depth can exist (occlusions, foreground objects, lighting effects etc) and must be taken care of. Our results are verified in a real television broadcasting studio with and without foreground dynamics. The proposed tracking solution is also compared with KinectFusion which is a recently presented alternative technique [START_REF] Newcombe | Kinectfusion: Real-time dense surface mapping and tracking[END_REF]. This work optimizes the previous work with near real-time simulation results [START_REF] Tykkälä | Direct iterative closest point for real-time visual odometry[END_REF] into real-time application use. Due to keyframes, the local drift will be negligible without depth map registration.

Related work

In this section, various image content based camera tracking approaches are discussed. Our focus is on affordable, real-time methods which can be robust in studio environments without extensive scene manipulation.

Sparse feature based SLAM

Visual Simultaneous Localization and Mapping (vS-LAM) methods aim at estimating camera poses and structure concurrently from a real-time image sequence. The same scene geometry may be visible only for short time when the camera is moving. Therefore, an optimal pose and structure configuration are commonly solved within a sliding window of n recent frames using bundle adjustment [START_REF] Laboratorium | für Informationstechnologie, Hannover, Voodoo camera tracker: A tool for the integration of virtual and real scenes[END_REF][START_REF] Konolige | FrameSLAM: from bundle adjustment to realtime visual mappping[END_REF]. The problem is formulated by a cost function which measures 2D projection discrepancy of the current pose and the structure configuration. The minimization requires an initial guess, which is iteratively improved until the cost decreases below a threshold. The relative pose between two views are often inferred from the essential matrix, which can be estimated from 2D point correspondencies. The structure can be initialized using direct depth measurements or by triangulation from 2D projection points. This pose estimation process suffers from some drawbacks. Feature point extraction and matching is an error-prone process and requires an outlier rejection mechanism such as RANSAC. Extraction of high quality features such as SIFT can be expensive for real-time systems and still mismatches can not be fully prevented. For example, repetitive patterns and homogeneous textureless regions easily produce mismatches. Typically the estimation suffers from time evolving drift, which is often corrected by a loop closure mechanism. In loop closure, the current viewpoint is identified using a previously stored map. The cumulated pose error can then be divided evenly along the loop for removing the drift. Loop closure mechanisms are problematic in real-time AR applications unless they remove drift nearly in every frame, because bigger corrections are often visually disturbing. PTAM system by Klein is able to do loop closure for each frame, but suffers from error-prone feature extraction and matching process [START_REF] Klein | Parallel tracking and mapping for small ar workspaces[END_REF].

Pose tracking using planar surfaces

By introducing a priori information of surrounding 3D geometry into the estimation, the camera tracking becomes easier to solve. Several algorithms have been developed to track a camera based on a set of known 3D planar surfaces [START_REF] Silveira | An efficient direct approach to visual SLAM[END_REF]. The planarity of the environment allows describing the optical flow analytically via homography mapping, which can be decomposed into three-dimensional rotation and translation components. Homography mapping is described as a 3 × 3 matrix H, which is traditionally solved from point correspondencies. After successful initialization, the pose increments can also be obtained by using homography-based warping function which minimizes appearance based cost function. Accurate camera tracking is possible using visible planar markers, but unfortunately they are visually disturbing during television broadcasting.

Direct methods with dense 3D structures

The direct methods use image alignment methods for pose estimation instead of depending on 2D feature points [START_REF] Tykkälä | Direct iterative closest point for real-time visual odometry[END_REF][START_REF] Comport | Accurate quadri-focal tracking for robust 3d visual odometry[END_REF][START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF]. By avoiding the error-prone feature extraction and matching process, precision and robustness can be increased. Direct methods minimize the photometrical error between a reference point cloud and the current RGB image. The drift can be further reduced by matching also the current depth map values with the reference points [START_REF] Tykkälä | Direct iterative closest point for real-time visual odometry[END_REF]. A joint SLAM cost function has been formulated, which integrates structure component into the same cost function [START_REF] Tykkälä | A dense structure model for image based stereo SLAM[END_REF]. Steinbrücker et al. have recently proposed minimization of a photometrical cost [START_REF] Steinbrücker | Real-time visual odometry from dense RGB-D images[END_REF], but do not take the benefit of Mestimator, which can increase robustness with neglible computational cost. The image warping function is used to model the appearance change due to camera motion. The warping requires a homography or a relatively accurate 3D estimate of the surrounding environment geometry. 3D point appearance generally depends on the lighting conditions and the point of view, but when observing subsequent frames, the points approximately maintain the same color. DTAM applies the photometrical cost function for pose estimation, but also refines 3D structure by minimizing photometrical reprojection error [START_REF] Newcombe | Dtam: Dense tracking and mapping in real-time[END_REF]. DTAM is designed to operate with monocular camera in real-time, but it has several drawbacks, such as a special initialization procedure, discretization of depth values, not being compatible with dynamic scenes, heavy hardware requirements and parameter tuning with the correct surface smoothness. A RGB-D sensor such as Microsoft Kinect is better suited for direct pose estimation because it produces a real-time feed of RGB-D measurements which can be photometrically aligned [START_REF] Audras | Real-time dense RGB-D localisation and mapping[END_REF]. KinectFusion formulates pose estimation as depth map alignment problem using traditional Iterative Closest Point (ICP) [START_REF] Newcombe | Kinectfusion: Real-time dense surface mapping and tracking[END_REF]. The system also aims at maintaining small drift by registering the current depth maps with a voxel grid based 3D model, which is filtered over time by weighted averaging.

Dense RGB-D tracking

Our previous work tracked camera with respect to the RGB-D measurements which were recently captured by a RGB-D sensor. In this work, the previous frame reference is replaced by fixed keyframes which are stored in advance to avoid drift. The pose estimation is defined as a direct image registration task between the current image I : R 2 ⇒ R and the nearest keyframe K = {P * , c * , T * }, 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. The 4 × 4 matrix T * defines the keyframe pose relative to the reference coordinate system. T * maps points from the reference to the keyframe coordinate system. * is used to denote fixed reference variables.

Our goal is to find the correct pose increment T(x) which minimizes scalar cost function 1 2 e(x) T e(x), where the residual is defined by

e = I w P * ; T(x) T -c * . (1)
T = T g (T * ) -1 is the base transform, which encodes the relative pose between the initial guess T g and the keyframe pose T * . w(P; T) is the warping function which transforms and projects 3D points by

w(P; T) = w(P; R t 0 1 ) = KD(N(RP + t), α), (2)
where N(p) = (p 1 /p 3 , p 2 /p 3 ) dehomogenizes a point, D(p, α) performs lens distortion using parameters α ∈ R 5 (constant) and K is a 3 × 3 intrinsic matrix of the camera (constant). We use standard Caltech distortion model

D(p, α) = p(1 + α 1 r 2 + α 2 r 4 + α 5 r 6 ) + dx 1 (3)
where r 2 = p 2 1 + p 2 2 and dx is the tangential distortion

dx = 2α 3 p 1 p 2 + α 4 (r 2 + 2p 2 1 ) α 3 (r 2 + 2p 2 2 ) + 2α 4 p 1 p 2 . (4)
T(x) ∈ SE(3) forms a Lie group by exponential mapping [START_REF] Ma | An invitation to 3-D vision: from images to geometric models[END_REF]:

T(x) = e A(x) , A(x) = [ω] × υ 0 0 , (5)
where x = (ω, υ) ∈ R 6 encodes relative 3D motion between two camera poses. ω is the rotation axis and ω is the rotation angle in radians, and υ is the velocity twist. [ ] × produces a skew-symmetric matrix which performs a cross product. The matrix exponential produces 3×3 rotation matrix R and 3×1 translation vector t embedded in a 4 × 4 matrix. This cost function form does not require any predefined temporal correspondencies. The warping function will produce the dense correspondency when the minimization has converged.

Minimization

The inverse compositional approach [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF] is adopted for efficient minimization of the cost function. The cost is reformulated as c * (x) = I * w P * ; e -A(x) , c w = I w P * ; T e = c * (x)c w , where reference colors in c * are now a function of the inverse motion increment (Figure 2). The current warped intensities c w are produced by resampling I using T.

The residual e is then minimized by resampling the reference image I * . The benefit can be observed from the form of the Jacobian

J i j = ∂c(0) ∂x = ∇I * (w(P i ; I)) ∂w(P i ; I) ∂x j , (6)
where

∇I(p) = [ ∂I(p) ∂x ∂I(p)
∂y ]. Now J does not depend on T anymore and it can be computed only once for each K * for better performance. 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. c * k (0) are bi-linearly interpolated only once using reference image I * . 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 (7)
where the derivative is zero when

J T Jx = -J T e(0) . (8)
Because J T J is a 6 × 6 positive definite matrix, the inversion can be efficiently calculated using the Cholesky method. The matrix multiplication associativity enables collecting the estimated increment into the base transform by T k e A( x) ⇒ T k+1 .

Cost function smoothness

The first-order Taylor approximation is required for efficient minimization and therefore the cost function must be assumed to be locally smooth. Unfortunately, the exact cost function is non-linear and non-continuous due to arbitrary surface and appearance functions and the perspective projection. Bi-linear interpolation is used for continuous sampling of the images and smoothness is assumed between the pixels. The motion is parameterized by the matrix exponential, whose derivatives at zero are simply constant matrices ∂T(0) ∂x j = ∂A(0) ∂x j . Strong photometrical variations at different resolutions are the main cue of camera motion. The cost smoothness increases with the number of 3D points because noise deviations average out. Larger convergence domains can be obtained by using downsampled images, provided that aliasing effects do not dominate. By experimentation, 80 × 60 images provide rough estimates which can be refined using higher resolution images. Discontinuities in the cost function may occur when switching between different resolutions. Maintaining a fixed structure throughout the minimization is important and thus we do not downsample the geometry even though low resolution intensity images are used.

Rejecting outliers using M-estimator

The warping function models the appearance change due to camera motion when 3D points are static with Lambertian reflectance. Some of the points are considered outliers, because this model is not sufficient to explain their appearance. For example, intensity changes due to occlusions, lighting effects and geometry dynamics have not been modeled. Outliers can be detected efficiently by comparing the intensity and depth values at the warped points and the reference points. Inliers always have small error whereas outliers may have any error value. Increased robustness is obtained by damping the high error values out from the estimation (Figure 3). Instead of using a fixed threshold, a certain adaptation level to image brighness variations is useful. One method to determine the damping weights is by the Tukey weighting function

u k = |e k | c * median(e a ) , (9)
w k c =        (1 -( u k b ) 2 ) 2 if |u k | <= b 0 if |u k | > b , (10)
e z = Z w(P; T) -e T 3 T P 1 , (11)
without adaptation, where Z : R 2 ⇒ R is the depth map function of the current RGB image and e T 3 = (0, 0, 1, 0) is used to obtain the current depth value. 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 objects. The depth weights are determined by

w z k = max 1 -e 2 z (k)/τ 2 , 0 2 . ( 12
)
The weighted Gauss-Newton step is obtained by rewriting equation 8 as follows

J T WJx = -J T We, ( 13
)
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 weights are quadratic and therefore square root is not evaluated in practice.

Selecting reference points

The points used in the motion estimation can be freely selected from a reference RGB-D measurement. The points which do not contribute to the residual vector through linearization are not useful. Thus, only the points, which are associated with the greatest absolute values of the Jacobian should be selected [START_REF] Meilland | A spherical robot-centered representation for urban navigation[END_REF]. The selected points will then have strong optical flow ∂w ∂x > θ w and/or strong image gradient ∂I ∂p > θ g . We use a simpler selection criteria and focus on the set of points P * s which merely have strong image gradient. This method fits GPU better, because it does not require sorting of the Jacobian elements.

P * s = { P k | ∇I * (w(P k ; I)) > θ g , P k ∈ P * }, ( 14
)
where θ g is the threshold which selects p percent of all points. p depends on the computational capacity available and the amount of image edges in the application context.

Point distortion modeling vs image undistortion

By modeling the distortion in the warping function, the input images do not need to be undistorted. This is beneficial because image undistortion requires additional resampling of the image. Resampling degrades raw RGB data, because bi-linear interpolation fails to predict exact colors precisely. Also the number of distortion operations for a set of points will be smaller than full image resolution even though warping has to be done many times per frame.

Real-time implementation

Our method was implemented in Ubuntu Linux environment using open software tools, Microsoft Kinect and a commodity PC laptop hardware on which the method runs real-time. The cost minimization requires special attention when aiming at an efficient real-time implementation. Because the computational phases benefit from parallel computing, we implemented the full method (algorithm 1) on a commodity GPU. The same minimization algorithm can be used with dense incremental tracking and keyframe tracking. The implementation scales into multiple threads/cores. Only Cholesky inversion and matrix exponential computation are executed on a single GPU thread. for all iterations j = {1 . . . n L } do Reduce linear system (sec. 4.3) 8:

Solving linear system for x (sec. 4.3) 9:

T ⇐ Te A( x) (sec. 4.4)

10:

end for 11: end for

Warping

The warping function is applied in parallel to the selected P k ∈ P * . The cost function (eq. 1) is evaluated using bilinear interpolation. Additionally the points are transformed into IR view for evaluating the nearest depth value (eq. 11). Unnecessary resampling errors are avoided by evaluating the depth values directly in the IR view. The depths are used to generate weighting W z . The warping is illustrated in Figure 4.

M-estimator

Tukey based weights require the median of the error distribution, which easily becomes a bottleneck of using M-estimator. Therefore, we use the histogram technique (eq. 15) to find an approximate median. The distribution of the residual is represented using a 64-bin histogram and n is set to half residual length. histogram64 routine fits this purpose, because it is fast enough to be executed in every iteration and is specifically designed for NVIDIA video cards [START_REF] Podlozhnyuk | Histogram calculation in CUDA[END_REF]. By experiment, 64 bins seems to provide sufficient adaptation to different error distribution profiles.

Linear system reduction

The reduction of the linear system means computing J T J and J T e. This phase compresses the linear equations from n-dimensional to six-dimensional. Reduction is required for each iteration and therefore should be implemented efficiently. Matrix multiplication is often computed in parallel by dividing k dot products into separate threads. In this case, we have only 36 dot products, but our video card can manage 1024 threads in parallel. To gain maximal efficiency, we parallelize the computation in the dot product direction instead. Each dot product has n elements, where n is usually 8192 or more. Dot products are reduced by dividing them into subblocks and summing them efficiently in parallel. The total sum is finally cumulated from the subsums. J T J is a positive semidefinite matrix which is also symmetric. This means it is sufficient to compute only the upper triangle of the values and mirror the results into the lower triangle. This property reduces the number of dot products from 36 to 21. J T e requires six dot products. The process is illustrated in Figure 5. The inversion can be efficiently calculated using Cholesky decomposition due to positive definite property. The explicit inverse of J T J can also be avoided by using the conjugate gradient method with six iterations. Both approaches are fast due to the small matrix size. The inversion is executed on a single GPU thread.

Evaluating matrix exponential

Matrix exponentials can be computed in closed form using the translation extended Rodriquez formula [START_REF] Ma | An invitation to 3-D vision: from images to geometric models[END_REF]. Unfortunately it is not numerically stable in the most important small angle case. Expokit, however, provides various numerical approaches which are guaranteed to produce smooth and precise mapping [START_REF] Sidje | Expokit: Software package for computing matrix exponentials[END_REF]. By comparing them with MATLAB default implementation, it seems that complex matrix exponential (zgpadm) is the most accurate. Thus, we choose to evaluate T(x) = e A(x) by the matrix exponential of a general complex matrix in full, using the irreducible rational Padé approximation combined with scaling-and-squaring. The exponential is computed in a single GPU thread due to sequential nature of the operation.

Selecting points on GPU

Even though keyframe reference points can be selected in a pre-process, a fast implementation is useful when executing dense tracking incrementally. We adopt the histogram technique to find the best points efficiently [START_REF] Tykkälä | Direct iterative closest point for real-time visual odometry[END_REF]. We seek such histogram bin B t for which

Σ 255 k=B t -1 h(k) < n <= Σ 255 k=B t h(k), ( 15
)
where n is the number of points to be selected and

G = |∇ u I * (w (P * ; I)) | + |∇ v I * (w (P * ; I)) | h = histogram(G/2.0)
The value range is bounded to [0, 255] which fits histogram256 method designed for NVIDIA video cards [START_REF] Podlozhnyuk | Histogram calculation in CUDA[END_REF]. The final n indices are collected into a packed array. Packing is important because it allows eliminating all non-interesting points from further processing in the pipeline. Packing is implemented on a GPU by assigning each thread a slice of the original index range to compress (Figure 6). Each thread stores the interesting points into a temporary packed buffer with the count. The counts must be collected from all threads to determine the final output range of each thread. This allows parallel point selection which scales into multiple threads. Due to the discretization into bins, n pixels may not match the bin boundary automatically. n should be a multiple of 1024 to match the maximum amount of threads on a GPU. This is why it is useful to classify points into

P semi = { |∇c k | = B t | P k ∈ P * } P good = { |∇c k | > B t | P k ∈ P * },
where P good are all selected and the remaining points are chosen from P semi to obtain exactly n selected points.

Preprocessing RGB images

The preprocessor converts 640×480 Microsoft Kinect Bayer images into a pyramid of 320 × 240, 160 × 120, and 80 × 60. A 5 × 5 Gaussian filter is used with downsampling of the high resolution images. Downsampling is almost lossless, because the Bayer images are redundant. 2 × 2 block averaging is used to produce the rest of the layers. The lower resolution x and y coordinates are obtained by

x L = 1 2 L x + 1 2 L+1 -1 2
, where L is the amount of layers in-between. The RGB pre-processing steps are per-pixel operations which are executed in separate threads on a GPU.

Microsoft Kinect Calibration

Microsoft Kinect is conveniently calibrated using a specialized toolbox by Herrera et al., which jointly estimates all calibration parameters [START_REF] Herrera | Joint depth and color camera calibration with distortion correction[END_REF]. A chessboard pattern is printed and a set of matching raw disparity images and RGB images are captured and loaded into the toolbox. The toolbox then semi-automatically provides the intrinsic matrices K IR , K RGB , the baseline transform T b between the views, reconstruction parameters (c 0 , c 1 ) ∈ R 2 and distortion coefficients α RGB , α 0 , α 1 and β. The raw disparity map D(u, v) is then undistorted by [START_REF] Podlozhnyuk | Histogram calculation in CUDA[END_REF] and converted into depth values by

u(d, α 0 , α 1 , β) = d + β(u, v) exp (α 0 -α 1 d),
z(d) = 1 c 0 + c 1 u(d) . (17)
The disparity undistortion model corrects exponential decay of accuracy, but also models per pixel distortions using a map β(u, v). The RGB image lens distortion is modeled using the standard Caltech parameters α RGB ∈ R 5 , which models radial distortions using three coefficients and tangential distortions using two coefficients [START_REF] Bouguet | Camera calibration toolbox for matlab[END_REF]. Typically only the first two radial components are necessary, and the rest can be fixed to zero to speed up computation. The toolbox allows fixing distortion parameters prior to calibration.

Generating point cloud from raw disparity map

The 3D points are generated in parallel by a baseline transform ). The intensity vector c * matching with the 3D points P k ∈ P is produced by bi-linear interpolation, because the points do not match with RGB image pixels.

P k = T b z(d k )K -1 IR (p k 1) T , (18)

Results

The minimization algorithm is implemented on a low-end NVIDIA NVS4200m GPU using CUDA. The GPU has 48 CUDA cores and it allows executing 1024 parallel threads on a single multi-streaming processor. Because our implementation divides the computational task into n threads, it is scalable and benefits from GPU hardware development. Currently the most powerful video cards, such as the NVIDIA Tesla K10, have 3072 CUDA cores. The computation time of the processing phases is illustrated in Figure 8. The green bars represent the minimization phases which directly scale into n threads and therefore become faster with more powerful GPU. 8192 points are selected and the minimization uses three multi-resolution layers 80×60, 160×120 and 320 × 240. The corresponding iteration counts are 2, 3, and 10. In effect, the system operates at 30Hz and the computation takes 26ms which leaves 7ms for rendering the augmented graphics.

Incremental Dense Tracking Experiments

The incremental dense tracking is sketched in Algorithm 2. In this section the dense tracking accuracy is evaluated 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]. Kinfu is the open source implementation of KinectFusion [START_REF] Rusu | 3D is here: Point Cloud Library (PCL)[END_REF] which is a recent technique for low-cost camera tracking and reconstruction. The proposed dense RGB-D tracking accuracy is compared with Kinfu (Figure 9). Kinfu experiments were executed on a powerful workstation GPU (NVIDIA Quadro 2000, 1024MB, 192

T re f = T cur 8:
end if 9: end for Table 5 shows the comparison between our method and Kinfu numerically using two Freiburg sequences with known ground truth trajectory. Our dense tracking drifts 1.08cm/s with the slower freiburg2/desk sequence and 2.60cm/s with the faster freiburg1/desk sequence. Kinfu has smaller drift with small voxel volumes (such as (3m) 3 ) when the scene geometry is versatile, but seems to suffer from gross tracking failures with bigger volumes such as (5m) 3 and (8m) 3 . Operation volume is limited because 512 3 voxel grid becomes too coarse and bigger grids not not fit into GPU memory. Also ICP breaks down easily when the scene contains mostly planar surface (e.g. floor) 1 . 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. One second average error was computed from the subsegment with median error to average out random perturbations and to neglect Kinfu tracking failures, which occurred in all cases except on freiburg2/desk using (3m) 3 volume. Our incremental dense tracking has generally smaller drift with the faster freiburg1/desk sequence, but both Kinfu (with the most compatible grid) and our method lost track once during the sequence. In this case (3m) 3 grid does not contain versatile geometry and the result is worse than with bigger volumes. Re-localization issues are discussed later in the paper. Our method has more delay than previously because the reference frames must be updated more frequently with the faster sequence. Kinfu's dependency on careful setting of volume size, and dependency on geometrical variations makes it unsuitable to be used in our application. In television production studios, scenes can easily be larger than (3m) 3 and sufficient geometrical variation is more difficult to guarantee than sufficient texturing. Our dense tracking operates without failures even when planar surfaces are present, because our cost function matches also texturing. Memory consumption can be low even in larger operating volumes, because the keyframe placement can be optimized based on the camera motion zones.

Dense RGB-D tracking using keyframes

Despite that our dense tracker works in cases where KinectFusion fails, the tracking suffers from timeevolving drift, which is not acceptable in studio use, because it causes virtual items to move away from their 1 Video: http://youtu.be/tNz1p1sdTrE Figure 9: Kinfu performance compared with incremental dense tracking using freiburg desk2 sequence with motion capture ground truth (green trajectory). Kinfu is executed with (3m) 3 , (5m) 3 and (8m) 3 voxel volumes. The red trajectories on the left are output from from Kinfu. Kinfu gains lower drift due to structure integration, but planar surfaces cause tracking failures. (3m) 3 volume does not contain floor and therefore Kinfu works well. On the right, the yellow trajectory is the proposed incremental RGB-D tracking result. Problems with planar surfaces do not exist, and the method allows larger operating volumes due to lower memory consumption. By using keyframes, the drift can be completely eliminated. The green dots reveal the selected points.

Algorithm 3 Keyframe tracking algorithm. T cur ⇐ TT * 5: end for correct pose [START_REF] Tykkälä | Direct iterative closest point for real-time visual odometry[END_REF]. A fixed 3D model is required to avoid drift (Figure 1). A keyframe model can be generated with various offline/online techniques prior to broadcasting [START_REF] Laboratorium | für Informationstechnologie, Hannover, Voodoo camera tracker: A tool for the integration of virtual and real scenes[END_REF][START_REF] Henry | RGB-D mapping: Using Kinect-style depth cameras for dense 3D modeling of indoor environments[END_REF]. The keyframe tracking is sketched in Algorithm 3. Tracking is initialized at the first keyframe whose T = I. The cost function (eq. 1) is then minimized in each frame to obtain 3D camera pose. The reference keyframe is switched when the current pose becomes closer to another keyframe. Keyframe tracking is driftless and very fast due to keyframe pre-computations and the utilization of GPU.

Keyframe model generation

In a studio environment, online structure correction mechanism is an unnecessary because the correct 3D model can be fixed prior to a broadcast. We are in-terested in rapid and semi-automatic reconstruction approaches, which allow verifying and editing the 3D model prior to broadcasting. We sweep the scene using RGB-D sensor and concurrently store the pose parameters and depth maps using our real-time dense tracking approach (Algorithm 2). At this time, the studio scene must not contain moving actors but merely static geometry. The final keyframes are incrementally chosen from the RGB-D stream by discarding the frames which are too close to the existing database (Section 6.2). A loopclosing mechanism could be used with more complex sweeps where the camera re-enters same view points. Loop closure can be implemented by alternating between Algorithms 2 and 3 and adding a keyframe only when necessary. Both methods are rapid and they produce a sub-optimal set of keyframes, which can be further refined in offline when necessary.

A simple keyframe model generation process is Because the keyframes are sparsely selected from the stream of RGB-D measurements, they do not utilize all information available. The intermediate point clouds, which are not selected as keyframes, are warped into a nearest keyframe and the final maps are filtered in postprocessing. In effect, this improves depth map accuracy and fills holes. A large portion of the outlier points can be neglected based on color deviation to the keyframe pixels. The initial depth map is obtained as the median of the warped depths similar to Hirschmüller [START_REF] Hirschmüller | Stereo processing by semi-global matching and mutual information[END_REF]. The inverse depth samples within z m ± δ are averaged to find the robust estimate (Figure 10). Inverse depths typically have Gaussian distribution in stereo based settings. δ is the depth sample window, which is user-specified. δ depends on the RGB-D sensor depth noise level. This procedure is comparable to the weighted sum of Truncated Signed Distance Functions [START_REF] Newcombe | Kinectfusion: Real-time dense surface mapping and tracking[END_REF], but the surface will always appear near the median distance, and therefore does not depend on the grid density. Sparse bundle adjustment can be used to remove the small global error at the end of estimated camera trajectory [START_REF] Lourakis | SBA: A Software Package for Generic Sparse Bundle Adjustment[END_REF]. SBA does not converge without a good initial guess. Sometimes an initial guess can be succesfully extracted from an unordered image set by extracting and matching feature points and initializing camera configuration using their geometrical relations. The process is error-prone because image content such as homogeneous texturing and repetitive patterns prevents reliable extraction/matching. The studio scenes in our case do not contain much texture variations, and therefore automatic feature extraction fails. In this case, however, dense tracking provides very good initial guess of the pose configuration. To be able to execute SBA, few 3D points and their 2D correspondence are required from the user. An easy-to-use interactive editor was implemented for this purpose. Multiple geometries disappear when executing SBA (Figure 11). The dense scene geometry becomes precise with correct poses. Notice especially the sofa shape before and after the correction.

Finding nearest keyframe

A similarity metric is required to find the nearest key pose T k . The challenge is to unify rotation and translation differences, because they are expressed in different units. First, we define the relative transformation

∆T k = T cur T -1 k ⇒ (θ k , v k , d k ), ( 19
)
where T cur is the current camera pose. The relative rotation is decomposed into angle-axis representation (θ, v) based on the Euler's rotation theorem. The translation between the poses is expressed as the vector d between the origins. We define a potential keyframe index set as

Ω = { |θ k | < θ max , ||d k || < d max || k ∈ [1, n] }, ( 20
)
where n is the number of keyframes in the database.

Thus Ω contains a subset of keyframe indices whose angle and distance are below user defined thresholds θ max

where K s is the nearest keyframe. In contrast to frustum intersection, this metric works also when the camera is rotating around z-axis. The test point set P is a sparse representation of the view frustum (Figure 12). In particular, the frustum is approximated by three sparse 2D grids, each having uniformly sampled depth coordinates in the overall depth range of the RGB-D sensor. The 3D points at different layers are generated by discrete steps along the viewing rays.

Dense keyframe tracking experiments

In Figure 13, it is shown how the drift increases with dense tracking when moving back and forth along a fixed rail in a studio environment. We measure distance to the ground truth in every frame as error metric. The ground truth is generated without bundle adjustment, simply by fitting a line based on the first sweep (alg. 2). Incremental tracking collects pose error from two different sources: 1) numerical inaccuracies in the pose estimation, 2) moving objects in the reference view. The drift problem is solved by tracking relative to keyframes which contain merely static background geometry (alg.

3). The demonstration video shows the difference between dense tracking with and without keyframes 3 .

3 http://youtu.be/zfKdZSkG4LU Figure 14 shows how the online tracking accuracy depends on the number of keyframes. The sequence is illustrated in Figure 1, but we measure keyframe switching error in an empty scene for eliminating occlusion effects. The ground truth is generated by applying bundle adjustment to 27 keyframes which are initialized by dense tracking. Keyframe tracking is then executed with a sparsified number of ground truth keyframes (14, 9, 7, 5) and the camera pose is compared in each frame to the corresponding ground truth pose. In long-term use, even a small number of keyframes eventually outperforms the dense tracking due to drift. However, a small keyframe number produces local drift which appears as error ramps at keyframe switching points. This can be visually disturbing. With sufficient number of keyframes, the error remains small. Keyframe pose consistency finally depends on the on the accuracy of bundle adjustment. In this experiment, we obtain larger keyframe switching effects due to bundle adjustment phase. Bundle adjustment removes global error but can introduce local variance to the keyposes if the 2D point correspondencies are not precise. With the demonstration video, 14 keyframes produces small switching effects.

Rejecting dynamic foreground in studio

During broadcasting, the scene will contain moving actors. The pose estimation should only use static background in order to avoid estimation bias. After the static keyframe model has been captured, the scene can be segmented into static and dynamic components on the current depth and intensity difference to the static keyframe. The weighting computed in section 3.3 enforces the estimation to use static background. The weighting is illustrated in Figure 15. When relying on incremental tracking (alg. 2), the reference RGB-D measurement may unfortunately contain also moving objects which interfere tracking. The effect can be noticed during the last cycles (Figure 13) where an actor is moving in the scene. The estimation requires sufficient amount of visible points and, therefore, foreground actors are not allowed to occlude more than a fraction of the selected points. Figure 13 illustrates how the moving actor does not disturb camera tracking when relying on keyframes.

Studio lighting conditions with Microsoft Kinect

In studio environments, the color constancy assumption works well, because lighting can be fixed and specular surfaces avoided. This implies that keyframe database is valid as long as the scene configuration is fixed. Interfering IR light which is not originated from the RGB-D sensor can easily be cut-off in studio environments. One problem, however, is that the Microsoft Kinect RGB camera is of low quality and in standard lighting conditions the colors can saturate (see Fig. 16). With professional HD cameras this is not the case since they have remarkably better dynamic range. To circumvent this limitation, the RGB-D sensor depth maps should be directly calibrated with a HD camera [START_REF] Herrera | Joint depth and color camera calibration with distortion correction[END_REF].

Constraints

If the camera is moving too quickly, the minimization may not converge, because a local optimization strategy is used. Global minimization strategies of the cost function are not discussed in this paper because they are computationally too expensive to operate in real-time. Motion capture systems, on the otherhand, are expensive but could be used to provide initial guess for a pose. To avoid re-localization needs completely, we operate with limited camera speed, sufficient scene texturing, sufficient keyframe density, and ensure that sufficient amount of selected keyframe points are visible despite the occluding actors.

Conclusions

In this work, an affordable real-time matchmoving solution has been developed, which can be used to produce broadcasts with interactive digital components, such as virtual characters and stage items. The solution performs in real-time using a RGB-D sensor, and a laptop with a low-end GPU. This system was able to accurately and robustly track camera in broadcasting studio sized operating volumes where the state-of-the-art approach fails. Time-evolving drift problem is solved by tracking the camera relative to the nearest keyframe. The keyframes were generated using incremental dense tracking, and fine-tuned using sparse bundle adjustment. M-estimator was enhanced by segmentationbased weights, which allows actors to move in the foreground while tracking the camera. When operating in the RGB-D sensor range, our pose estimation accuracy depends mostly on texturing, which is trivial to manipulate/add in studio environments. Therefore the proposed method is easily adaptable to different types of studio scene settings, whereas KinectFusion is constrained by the amount of geometrical variations in the scene. Camera tracking has been demonstrated in a real broadcast studio with and without dynamic components. Drift and keyframe switching errors have also been quantified. Future work will address the practical issues of how studio staff and cameramen can use this computer vision system in live broadcasts. We will calibrate RGB-D sensor with a professional HD camera for better image quality. Moreover, combination of the best properties of our approach and KinectFusion will be investigated.

Figure 2 :

 2 Figure 2: Inverse compositional image alignment. The inverse increment is found which minimizes the difference between I w and I * . Linearization is specially computed only once using I * and P * .

  where e a = (|e 1 |, |e 2 |, . . . , |e n |) is a vector of absolute residual values, c = 1.4826 is the robust standard deviation, and b = 4.6851 is the Tukey specific constant. Thus, Tukey weighting generates adaptive weighting based on the statistical distribution of the residual.The depth correlation weights are computed directly from the depth residual

Figure 3 :

 3 Figure 3: Tukey M-estimator effect illustrated. The points with high error values are given small or zero weight, because their appearance change is not explained by the motion model.

Algorithm 1

 1 Minimization on GPU.Input: I L={1,2,3} with 320 × 240, 160 × 120, and 80 × 60 sizes.Z in 320 × 240. Iteration counts {n 1 , n 2 , n 3 }. T = T 0 . Output: Relative pose T.1: for all multi-resolution layers L = {3, 2, 1} do 2:

3 :

 3 Compute residual e and W z (sec. 4.1) 4: Determine M-estimator weights W c (sec. 4.2) 5: W ⇐ W c * W z

Figure 4 :

 4 Figure 4: The warping of points from keyframe into Microsoft Kinect RGB and IR image coordinate systems. The additional depth lookup from an IR image is required by Equations 1 and 11.

Figure 5 :

 5 Figure 5: Efficient reduction on GPU. Dot products are divided into sub-ranges which are reduced in parallel. J T J is symmetric 6 × 6 matrix and elements can be mirrored with respect to the diagonal.

Figure 6 :

 6 Figure 6: Index buffer packing illustrated. Each thread selects points from a sub-range. The sub-counts are shared between the threads and the output buffer written.

Figure 7 :

 7 Figure 7: RGB and disparity images, and distortion pattern β [18]

Figure 8 :Algorithm 2

 82 Figure 8: GPU computation time for one frame. Green bars represent phases which are scalable to n threads.The red bar must be executed on a single thread. The gray bar represents the CUDA interop delay of buffer locking which is required for visualization only. d2z is the disparity map to zmap conversion, layers is the conversion to grayscale multi-resolution format, pcloud is the 3D reconstruction, warp contains the cost function evaluations, mest is the M-estimation, reduct is the reduction of the linear system, and minc contains the delay of solving the motion parameters and evaluating the matrix exponential.

Require:

  Keyframe database available Input: T cur = I Output: Trajectory {I, T 1 cur , . . . , T n cur }. 1: for each RGB-D measurement do 2: {P * , c * , T * } ⇐FindKeyframe(T cur ) 3:T ⇐Minimize(T cur (T * ) -1 , P * , c * , I cur ) 4:

•

  Incremental dense tracking (short sweep) Followed by the optional offline refinement phases • Depth map fusion • Sparse bundle adjustment (semi-automatic annotation of keypoints)

Figure 10 :

 10 Figure 10: Robust depth estimation near median depth illustrated.

Figure 11 :

 11 Figure 11: The effect of sparse bundle adjustment illustrated. Minor error at the end of the trajectory is corrected in the right side image. The dense scene geometry becomes precise with correct poses. Notice especially the sofa shape before and after the correction.

Figure 12 :

 12 Figure 12: The 3D test point set P used to compare camera poses is illustrated. P approximates the view frustum by three 2D point grids.The projection error of P between the current view and each keyframe view is compared to find the nearest keyframe.

Figure 13 :

 13 Figure13: Camera moving forward and backward along a fixed 3.30m studio rail. The top figure shows three images taken from the beginning, middle and the end of the rail. Green region represents the selected points. On the bottom, the comparison between incremental tracking and keyframe tracking. In incremental tracking drift increases in time, but keyframe tracking maintains small bounded error. A person is moving in the scene during the last cycles.

Figure 14 :

 14 Figure 14: The translation and rotation error as a function of keyframe count. The sequence is illustrated in Figure 1. Keyframe switching error depends on the amount of keyframes. Over 14 keyframes are required to obtain sufficiently small keyframe switching error.

Figure 15 :

 15 Figure 15: The points which currently participate in pose estimation are selected based on the image gradient magnitude. The dynamic foreground is given zero weighting in the estimation based on large depth and intensity difference to the static keyframe. The weight 1.0 is assigned to the green points and 0.0 to the red points.

Figure 16 :

 16 Figure 16: HD camera (left) and Microsoft Kinect RGB image (right) in typical studio lighting conditions.

Table 1 :

 1 Drift and delay compared to Kinfu with (3m) 3 and (8m)3 operating volume.

	Dataset	Our drift	Kinfu(3)	Kinfu(8)	Motion
	freiburg1/desk	2.60cm/s	8.40cm/s	3.97cm/s 41.3cm/s
		52.2ms	135ms	135ms
	freiburg2/desk	1.08cm/s	0.64cm/s	1.30cm/s 19.3cm/s
		35.5ms	135ms	135ms

Acknowledgements

We thank Heikki Ortamo and Jori Pölkki for their professional support in a TV brodcasting studio.


O Golinelli

K Mallick

Connected Operators for the Totally Asymmetric Exclusion Process

Keywords: Non-Equilibrium Statistical Mechanics, ASEP, Exact Results, Algebraic Bethe Ansatz Pacs numbers: 02.30.Ik, 02.50.-r, 75.10.Pq

We fully elucidate the structure of the hierarchy of the connected operators that commute with the Markov matrix of the Totally Asymmetric Exclusion Process (TASEP). We prove for the connected operators a combinatorial formula that was conjectured in a previous work. Our derivation is purely algebraic and relies on the algebra generated by the local jump operators involved in the TASEP.

Introduction

The Asymmetric Simple Exclusion Process (ASEP) is a lattice model of particles with hard core interactions. Due to its simplicity, the ASEP appears as a minimal model in many different contexts such as one-dimensional transport phenomena, molecular motors and traffic models. From a theoretical point of view, this model has become a paradigm in the field of non-equilibrium statistical mechanics; many exact results have been derived using various methods, such as continuous limits, Bethe Ansatz and matrix Ansatz (for reviews, see e.g., Spohn 1991, Derrida 1998, Schütz 2001, Golinelli and Mallick 2006).

In a recent work (Golinelli and Mallick 2007), we applied the algebraic Bethe Ansatz technique to the Totally Asymmetric Exclusion Process (TASEP).

This method allowed us to construct a hierarchy of 'generalized Hamiltonians' that contain the Markov matrix and commute with each other. Using the algebraic relations satisfied by the local jump operators, we derived explicit formulae for the transfer matrix and the generalized Hamiltonians, generated from the transfer matrix. We showed that the transfer matrix can be interpreted as the generator of a discrete time Markov process and we described the actions of the generalized Hamiltonians. These actions are non-local because they involve non-connected bonds of the lattice. However, connected operators are generated by taking the logarithm of the transfer matrix. We conjectured for the connected operators a combinatorial formula that was verified for the first ten connected operators by using a symbolic calculation program.

The aim of the present work is to present an analytical calculation of the connected operators and to prove the formula that was proposed in (Golinelli and Mallick 2007). This paper is a sequel of our previous work, however, in section 2, we briefly review the main definitions and results already obtained so that this work can be read in a fairly self-contained manner. In section 3, we derive the general expression of the connected operators.

Review of known results

We first recall the dynamical rules that define the TASEP with n particles on a periodic 1-d ring with L sites labelled i = 1, . . . , L. The particles move according to the following dynamics: during the time interval [t, t + dt], a particle on a site i jumps with probability dt to the neighboring site i + 1, if this site is empty. This exclusion rule which forbids to have more than one particle per site, mimics a hard-core interaction between particles. Because the particles can jump only in one direction this process is called totally asymmetric. The total number n of particles is conserved. The TASEP being a continuous-time Markov process, its dynamics is entirely encoded in a 2 L × 2 L Markov matrix M, that describes the evolution of the probability distribution of the system at time t. The Markov matrix can be written as

M = L i=1 M i , (1)
where the local jump operator M i affects only the sites i and i + 1 and represents the contribution to the dynamics of jumps from the site i to i + 1.

The TASEP algebra

The local jump operators satisfy a set of algebraic equations :

M 2 i = -M i , (2) M i M i+1 M i = M i+1 M i M i+1 = 0, (3) [M i , M j ] = 0 if |i -j| > 1.
(4)

These relations can be obtained as a limiting form of the Temperley-Lieb algebra. On the ring we have periodic boundary conditions :

M i+L = M i .
The local jumps matrices define an algebra. Any product of the M i 's will be called a word. The length of a given word is the minimal number of operators M i required to write it. A word, that can not be simplified further by using the algebraic rules above, will be called a reduced word. Consider any word W and call I(W ) the set of indices i of the operators M i that compose it (indices are enumerated without repetitions). We remark that, if W is not annihilated by application of rule (3), the simplification rules (2, 4) do not alter the set I(W ), i.e., these rules do not introduce any new index or suppress any existing index in I(W ). This crucial property is not valid for the algebra associated with the partially asymmetric exclusion process (see Golinelli and Mallick 2006).

Using the relation (2) we observe that for any i and any real number λ = 1 we have

(1 + λM i ) -1 = (1 + αM i ) with α = λ λ -1 . (5)

Simple words

A simple word of length k is defined as a word M σ(1) M σ(2) . . . M σ(k) , where σ is a permutation on the set {1, 2, . . . , k}. The commutation rule (4) implies that only the relative position of M i with respect to M i±1 matters. A simple word of length k can therefore be written as W k (s 2 , s 3 , . . . , s k ) where the boolean variable s j for 2 ≤ j ≤ k is defined as follows :

s j = 0 if M j is on the left of M j-1 and s j = 1 if M j is on the right of M j-1 . Equivalently, W k (s 2 , s 3 , . . . , s k ) is uniquely defined by the recursion relation W k (s 2 , s 3 , . . . , s k-1 , 1) = W k-1 (s 2 , s 3 , . . . , s k-1 ) M k , (6) W k (s 2 , s 3 , . . . , s k-1 , 0) = M k W k-1 (s 2 , s 3 , . . . , s k-1 ) . ( 7
)
The set of the 2 k-1 simple words of length k will be called W k . For a simple word W k , we define u(W k ) to be the number of inversions in W k , i.e., the number of times that M j is on the left of M j-1 :

u(W k (s 2 , s 3 , . . . , s k )) = k j=2 (1 -s j ) . (8)
We remark that simple words are connected, they cannot be factorized in two (or more) commuting words.

Ring-ordered product

Because of the periodic boundary conditions, products of local jump operators must be ordered adequately. In the following we shall need to use a ring ordered product O () which acts on words of the type

W = M i 1 M i 2 . . . M i k with 1 ≤ i 1 < i 2 < . . . < i k ≤ L , (9)
by changing the positions of matrices that appear in W according to the following rules :

(i) If i 1 > 1 or i k < L, we define O (W ) = W . The word W is well- ordered.
(ii) If i 1 = 1 and i k = L, we first write W as a product of two blocks, W = AB, such that B = M b M b+1 . . . M L is the maximal block of matrices with consecutive indices that contains M L , and A = M 1 M i 2 . . . M ia , with i a < b -1, contains the remaining terms. We then define

O (W ) = O (AB) = BA = M b M b+1 . . . M L M 1 M i 2 . . . M ia . (10)
(iii) The previous definition makes sense only for k < L. Indeed, when k = L, we have W = M 1 M 2 . . . M L and it is not possible to split W in two different blocks A and B. For this special case, we define

O (M 1 M 2 . . . M L ) = |1, 1, . . . , 1 1, 1, . . . , 1| , (11)
which is the projector on the 'full' configuration with all sites occupied. The ring-ordering O () is extended by linearity to the vector space spanned by words of the type described above.

Transfer matrix and generalized Hamiltonians H k

The algebraic Bethe Ansatz allows to construct a one parameter commuting family of transfer matrices, t(λ), that contains the translation operator T = t(1) and the Markov matrix M = t ′ (0). For 0 ≤ λ ≤ 1, the operator t(λ) can be interpreted as a discrete time process with non-local jumps : a hole located on the right of a cluster of p particles can jump a distance k in the backward direction, with probability λ k (1 -λ) for 1 ≤ k < p, and with probability λ p for k = p. The probability that this hole does not jump at all is 1 -λ. This model is equivalent to the 3-D anisotropic percolation model of Rajesh and Dhar (1998) and to a 2-D five-vertex model. It is also an adaptation on a periodic lattice of the ASEP with a backwardordered sequential update (Rajewsky et al. 1996, Brankov et al. 2004), and equivalently of an asymmetric fragmentation process (Rákos and Schütz 2005).

The operator t(λ) is a polynomial in λ of degree L given by

t(λ) = 1 + L k=1 λ k H k , (12)
where the generalized Hamiltonians H k are non-local operators that act on the configuration space. [We emphasize that the notation used here is different from that of our previous work : t(λ) was denoted by t g (λ) in (Golinelli and Mallick 2007).]

We have H 1 = M and more generally, as shown in (Golinelli and Mallick 2007), H k is a homogeneous sum of words of length k

H k = 1≤i 1 <i 2 <...<i k ≤L O (M i 1 M i 2 . . . M i k ) , (13)
where O () represents the ring ordered product that embodies the periodicity and the translation-invariance constraints.

For a system of size L with N particles only H 1 , H 2 , . . . , H N have a nontrivial action. Because we are interested only in the case N ≤ L -1 (the full system as no dynamics) there are at most L -1 operators H k that have a non-trivial action.

3 The connected operators F k

Definition

The generalized Hamiltonians H k and the transfer matrix t(λ) have non-local actions and couple particles with arbitrary distances between them. Besides H k is a highly non-extensive quantity as it involves generically a number of terms of order L k . As usual, the local connected and extensive operators are obtained by taking the logarithm of the transfer matrix. For k ≥ 1, the connected Hamiltonians F k are defined as

ln t(λ) = ∞ k=1 λ k k F k . ( 14
)
Taking the derivative of this equation with respect to λ and recalling that t(λ) commutes with t ′ (λ), we obtain

∞ k=1 λ k F k = λ t(λ) -1 t ′ (λ) . ( 15
)
Expanding t(λ) -1 with respect to λ, this formula allows to calculate F k as a polynomial function of H 1 , . . . , H k . For example Golinelli and Mallick 2007). By using (13), we observe that F k is a priori a linear combination of products of k local operators M i . However this expression can be simplified by using the algebraic rules (2, 3, 4) and in fine, F k will be a linear combination of reduced words of length j ≤ k.

F 1 = H 1 , F 2 = 2H 2 -H 2 1 , etc... (see
Because of the ring-ordered product that appears in the expression (13) of the H k 's, it is difficult to derive an expression of F k in terms of the local jump operators. An exact formula for the F k with k ≤ 10 was obtained in (Golinelli and Mallick 2007) by using a computer program and a general expression was conjectured for all k. In the following, the conjectured formula is derived and proved rigorously.

Elimination of the ring-ordered product

The expression λ k F k can be written as a linear combination of reduced words W . We know from formula (13) that at most L -1 operators H k are independent in a system of size L, we shall therefore calculate F k only for k ≤ L -1. Thus, we need to consider reduced words of length j ≤ L -1. Let W be such a word, and I(W ) be the set of indices of the operators M i that compose W ; our aim is to find the expression of W and to calculate its prefactor from equation ( 15). Because the rules (2, 4) do not suppress or add any new index, the following property is true : if a word W ′ appearing in λ t(λ) -1 t ′ (λ) is such that I(W ′ ) = I(W ) then even after simplification, W ′ will remain different from W . Therefore, the prefactor of W in λ k F k is the same as the prefactor of W in

λ t I (λ) -1 t ′ I (λ) where t I (λ) = O i∈I (1 + λM i ) with I(W ) ⊂ I . ( 16
)
Because F k commutes with the translation operator T , then for any r = 1, . . . , L-1, the prefactor of W = M i 1 M i 2 . . . M i j is the same as the prefactor of T r MT -r = M r+i 1 M r+i 2 . . . M r+i j . Furthermore, any word W of size k ≤ L -1 is equivalent, by a translation, to a word that contains M 1 and not M L : indeed, there exists at least one index i 0 such that i 0 / ∈ I(W ) and (i 0 + 1) ∈ I(W ) and it is thus sufficient to translate W by r = L -i 0 .

In conclusion, it is enough to study in expression (15), the reduced words W with set of indices included in

I * = {1, 2, . . . , L -1} . ( 17
)
Because the index L does not appear in I * , the ring-ordered product has a trivial action in equation ( 16) and we have

t I * (λ) = (1 + λM 1 )(1 + λM 2 ) . . . (1 + λM L-1 ) . (18)
We have thus been able to eliminate the ring-ordered product.

Explicit formula for the connected operators

In equation ( 18), differentiating t I * (λ) with respect to λ, we have

t ′ I * (λ) = L-1 i=1 (1 + λM 1 ) . . . (1 + λM i-1 )M i (1 + λM i+1 ) . . . (1 + λM L-1 ) . (19)
Using equation ( 5) we obtain

t I * (λ) -1 = (1 + αM L-1 )(1 + αM L-2 ) . . . (1 + αM 1 ) , with α = λ λ -1 . (20)
Noticing that λ(1

+ αM i )M i = -αM i , we deduce λ t I * (λ) -1 t ′ I * (λ) = (21) -α L-1 i=1 (1 + αM L-1 ) . . . (1 + αM i+1 )M i (1 + λM i+1 ) . . . (1 + λM L-1 ) .
The ith term in this sum contains words with indices between i and L -1.

Because we are looking for the words that contain the operator M 1 , we must consider only the first term in this sum, which we note by

Q Q = -α(1 + αM L-1 ) . . . (1 + αM 2 )M 1 (1 + λM 2 ) . . . (1 + λM L-1 ) . ( 22
)
In the appendix, we show that

Q = R 1 + R 2 + . . . + R L-1 , (23)
where R i is defined by the recursion :

R 1 = -αM 1 , (24) R i = λR i-1 M i + αM i R i-1 for i ≥ 2 . ( 25
)
To summarize, all the words in ∞ k=1 λ k F k that contain M 1 and not M L are given by Q = R 1 +R 2 +. . .+R L-1 . From the recursion relation (25) we deduce that R i is a linear combination of the 2 i-1 simple words W i (s 2 , s 3 , . . . , s i ) defined in section 2.1. Furthermore, we observe from (25) that a factor λ appears if s i = 1 and a factor α = λ/(λ -1) appears if s i = 0. Therefore, the coefficient f (W ) of W = W i (s 2 , s 3 , . . . , s i ) in Q is given by

f (W ) = (-1) u λ i (1 -λ) u+1 = (-1) u ∞ j=0 u+j j λ i+j ( 26
)
where i is the length of W and u = u(W ) is its inversion number, defined in equation ( 8). We have thus shown that

Q = L-1 i=1 W ∈W i f (W ) W = L-1 i=1 W ∈W i W ∞ j=0 (-1) u(W ) u(W )+j j λ i+j , (27)
where W i is the set of simple words of length i.

Finally, we recall that the coefficient in ∞ k=1 λ k F k of a reduced word W that contains M 1 and not M L is the same as its coefficient in Q. Extracting the term of order λ k in equation ( 27) we deduce that any word W in F k that contains M 1 and not M L is a simple word of length i ≤ k and its prefactor is given by (-1) u(W ) u(W )+k-i k-i . The full expression of F k is obtained by applying the translation operator to the expression (27); indeed any word in F k can be uniquely obtained by translating a simple word in F k that contains M 1 and not M L . We conclude that for k < L,

F k = T k i=1 W ∈W i (-1) u(W ) k-i+u(W ) k-i W , ( 28
)
where T is the translation-symmetrizator that acts on any operator A as follows :

T A = L-1 i=0 T i A T -i .
The presence of T in equation ( 28) insures that F k is invariant by translation on the periodic system of size L. All simple words being connected, we finally remark that formula (28) implies that F k is a connected operator.

Conclusion

By using the algebraic properties of the TASEP algebra (2-4), we have derived an exact combinatorial expression for the family of connected operators that commute with the Markov matrix. This calculation allows to fully elucidate the hierarchical structure obtained from the Algebraic Bethe Ansatz. It would be of a great interest to extend our result to the partially asymmetric exclusion process (PASEP), in which a particle can make forward jumps with probability p and backward jumps with probability q. In particular, we recall that the symmetric exclusion process is equivalent to the Heisenberg spin chain : in this case the connected operators have been calculated only for the lowest orders (Fabricius et al., 1990). This is a challenging and difficult problem. In our derivation we used a fundamental property of the TASEP algebra : the rules (2-4) when applied to a word W either cancel W or conserve the set of indices I(W ). The algebra associated with PASEP violates this crucial property because there we have M i M i+1 M i = pq M i . Therefore the method followed here does not have a straightforward extension to the PASEP case.

Appendix: Proof of equation (23)

Let us define the following series

Q 1 = -αM 1 , ( 29
) Q i = (1 + αM i )Q i-1 (1 + λM i ) for i ≥ 2 . ( 30
)
We remark that Q defined in equation ( 22) is given by Q = Q L-1 . Let us consider R i defined by the recursion (25). The indices that appear in the words of Q i and R i belong to {1, 2, . . . , i}. Therefore, we have [R j , M i ] = 0 for j ≤ i -2 , (31)

because the operators M 1 , M 2 , . . . , M j that compose R j commute with M i . From equations ( 31) and ( 5), we obtain

(1 + αM i )R j (1 + λM i ) = R j for j ≤ i -2 . ( 32
)
Furthermore, from (25), we obtain

M i R i-1 M i = λM i R i-2 M i-1 M i + αM i M i-1 R i-2 M i . ( 33
)
Because M i commutes with R i-2 , we can use the relation

M i M i-1 M i = 0 to deduce that M i R i-1 M i = 0 . (34)

Using equation (34), we find

From equations ( 32) and ( 35), we prove that the (unique) solution of the recursion relation ( 30) is given by equation ( 23), • Spohn H., 1991, Large scale dynamics of interacting particles, Springer, New-York.