Vincent Acary

Guillaume James

Franck Pérignon

Periodic motions of coupled impact oscillators

We study existence of time-periodic oscillations in a chain of coupled impact oscillators, for rigid impacts without energy dissipation. We formulate the search of periodic solutions as a boundary value problem incorporating unilateral constraints. This problem is solved numerically and different solution branches corresponding to nonlinear localized modes (breathers) and normal modes are computed.

Understanding the dynamics of nonlinear lattices (i.e. large networks of coupled nonlinear oscillators) is a problem of fundamental importance in mechanics, condensed matter physics or biology. One of the major issues concerns the mathematical analysis and numerical computation of special classes of nonlinear time-periodic oscillations that organize the dynamics in many situations. In particular, spatially periodic waves (standing waves or periodic traveling waves) and spatially localized waves (breathers) are the object of intensive research [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF][START_REF] Flach | Discrete breathers : advances in theory and applications[END_REF]. In this context, many theoretical and numerical works have focused on smooth nonlinear systems, whereas relatively few mathematical results are available for waves in nonsmooth infinite lattices [START_REF] Gendelman | Discrete breathers in vibroimpact chains: analytic solutions[END_REF][START_REF] Gendelman | Exact solutions for discrete breathers in a forced-damped chain[END_REF][START_REF] Lebellego | Phénomènes ondulatoires dans un modèle discret de faille sismique[END_REF]. Developing theoretical and numerical tools for the analysis of nonlinear waves in nonsmooth systems is extremely important for applications, in particular in the context of impact mechanics where unilateral contacts and friction come into play. Spatially discrete lattice models are frequently encountered in this context, in particular for the modeling of waves in multibody mechanical systems (e.g. granular media) or in finite element models of continuum systems. A classical example illustrating the latter case concerns thin oscillating mechanical structures (a string under tension or a clamped beam) contacting rigid obstacles [START_REF] Schatzman | A Hyperbolic problem of second order with unilateral constraints : the vibrating string with a concave obstacle[END_REF][START_REF] Ahn | An Euler-Bernouilli beam with dynamic frictionless contact : penalty approximation and existence[END_REF]. Such a structure can be described by a one-dimensional finite-element model involving a large number of degrees of freedom. The contact force between the string/beam and a rigid obstacle is either measure-valued (for rebounds with velocity jumps at contact times) or set-valued (if a wrapping of the string on the obstacle occurs). Although nonlinear periodic waves are observed in experiments [START_REF] Astashev | Experimental investigation of vibrations of strings interacting with point obstacles[END_REF], relatively little is known from a theoretical point of view on their existence and stability. Existence results have been derived in particular cases, for a continuum string model with point-mass or plane obstacle ( [START_REF] Cabannes | Presentation of software for movies of vibrating strings with obstacles[END_REF] and references therein). In addition, the existence and stability of time-periodic breathers (spatially localized oscillations) has been analyzed for discrete linear chains with a single node undergoing rigid impacts, both for conservative systems [START_REF] Gendelman | Discrete breathers in vibroimpact chains: analytic solutions[END_REF] and forced systems with dissipative impacts [START_REF] Gendelman | Exact solutions for discrete breathers in a forced-damped chain[END_REF].

In this work, we study the existence and stability of time-periodic oscillations in a chain of linearly coupled impact oscillators reminiscent of a model analyzed in [START_REF] Gendelman | Discrete breathers in vibroimpact chains: analytic solutions[END_REF], for rigid impacts without energy dissipation. We introduce a numerical method allowing to compute branches of time-periodic solutions when an arbitrary number of nodes undergo rigid impacts. For this purpose, we reformulate the search of periodic solutions as a boundary value problem incorporating unilateral constraints. We illustrate this numerical approach by computing some families of nonlinear spatially localized modes (breathers) and extended modes. We consider an infinite chain of impact ocillators with positions described by an infinite vector y(t) ∈ ℓ ∞ (Z) (the space of bounded sequences on Z). The dynamics is described by the following complementarity system

ÿn + y n -γ (∆y) n = λ n , n ∈ Z, (1)
0 ≤ λ ⊥ (y + ½) ≥ 0, (2)
if ẏn (t -) < 0 and y n (t) = -1 then ẏn (t

+ ) = -ẏn (t -), (3)
where (∆y) n = y n+1 -2 y n + y n-1 defines a discrete Laplacian operator, ½ denotes the constant sequence with all terms equal to unity and γ ≥ 0 is a parameter. Non-dissipative impacts occur for y(t) = -1 and give rise to impulsive reaction forces λ(t). We look for T -periodic solutions even in time, and assume each particle undergoes at most one impact during each period of oscillation. Consequently, for a given particle, impacts either occur at half-period multiples or do not occur all. We denote by I k ⊂ Z with k = 1 or 2 the index sets of particles impacting at t = (2p + k) T /2 for all p ∈ Z, and by I 0 := Z \ (I 1 ∪ I 2 ) the index set corresponding to non-impacting particles. We have thus λ n = 0 for all n ∈ I 0 and

λ n = 2 ẏn ( kT + 2 ) p∈Z δ (p+ k 2 )T for all n ∈ I k .
Introducing the splitting y = (y (0) , y (1) , y (2) ) corresponding to Z = I 0 ∪ I 1 ∪ I 2 , the above system can be reformulated as a boundary value problem on a half-period interval (0, T /2),

ÿn + y n -γ (∆y) n = 0, n ∈ Z, t ∈ (0, T /2), (4)
with boundary conditions

ẏ(i) (0) = 0 for i ∈ I 0 ∪ I 1 , y (2) (0) = -½, ẏ(i) (T /2) = 0 for i ∈ I 0 ∪ I 2 , y (1) (T /2) = -½, (5)
and constraint

y(t) + ½ > 0, t ∈ (0, T /2). (6)
We solve this problem numerically for a chain of N = 100 oscillators with periodic boundary conditions. We use a shooting method, i.e. determine z = (y (0) (0), y (1) (0), ẏ(2) (0)) ∈ R N such that the three boundary conditions of ( 5) at t = T /2 are satisfied. This requires to solve a linear system for z obtained through time-integration of the linear ODE (4) (the case of nonlinear local or interaction potentials could be addressed similarly using a Newton method). The constraint ( 6) is checked a posteriori. Solution branches are continued for fixed values of T , varying the linear stiffness γ and starting from the uncoupled (or "anticontinuum") limit γ = 0 [START_REF] Flach | Discrete breathers : advances in theory and applications[END_REF]. In this limit, for all fixed T ∈ (π, 2π), a choice of impacting particles and phases (determined by I 1 , I 2 ) selects a unique solution which can be continued up to some maximal value of γ. The linear stability of periodic solutions is analyzed through the eigenvalues of an associated monodromy matrix.

To perform this computation, we integrate ( 1)-( 2)-(3) numerically using the Siconos software for nonsmooth dynamical systems [10]. As an example, we describe in figure 1 three families of periodic solutions obtained with this method, namely two breather solutions (site-centered or bond-centered) and a spatially extended solution (nonlinear normal mode). The computation of periodic solutions based on the above approach is much more effective than numerical continuation of periodic solutions based on compliant models. In the latter case, impacts are described by smooth nonlinear Hertzian type potentials leading to stiff ODE and costly numerical continuation. Future extensions of this work will include an analytical continuation and stability analysis based on the same approach, the inclusion of dissipative impacts and forcing, and the application of the method to more complex finite-element models of continuous systems under impacts.

Figure 1 :

 1 Figure1: Computation of different periodic solutions for T ≈ 4.7. The left column displays particle positions at t = 0 for γ = 0.16, for two breather solutions with I2 = {50}, I1 = ∅ (site-centered breather, top panel) and I2 = {49}, I1 = {50} (bond-centered breather, middle panel), and for a nonlinear normal mode with spatial period two and I0 = 2Z, I1 = ∅ (bottom panel). These solutions can be continued for γ ∈ [0, γmax) with γmax ≈ 0.19. The site-centered breather is linearly stable for γ < γc ≈ 0.13, after which it becomes unstable (the top right panel displays the moduli of the corresponding Floquet eigenvalues). The time evolution of the position and velocity of the impacting particle (n = 50) is illustrated over a few periods for γ = 0.16 (right column, middle and bottom panels). The bond-centered breather and nonlinear normal mode are both unstable for all values of γ.


Miles Hansard

Georgios Evangelidis

Quentin Pelorson

Radu Horaud
email: radu.horaud@inria.fr

Cross-calibration of Time-of-flight and Colour Cameras

Keywords: Camera networks, time-of-flight cameras, depth cameras, camera calibration, 3D reconstruction, RGB-D data

Introduction

The segmentation of multi-view video data, with respect to physically distinct objects of interest, is an essential task in automatic scene-interpretation. Visual segmentation can be based on colour, texture, parallax and motion information (e.g. [START_REF] Hoiem | Recovering surface layout from an image[END_REF][START_REF] Tron | A benchmark for the comparison of 3-d motion segmentation algorithms[END_REF]). The task remains very difficult, however, owing to the combined effects of non-rigid surfaces, variable lighting, and occlusion. It has become clear that depth cameras can make an important contribution to scene understanding, by enabling direct depth segmentation, based on the measured scene-structure (see CVIU special issue [START_REF] Savarese | Special issue on 3D representation for object and scene recognition[END_REF]). This approach is also highly effective for dynamic tasks, such as body tracking and action recognition [START_REF] Liu | Special issue on visual understanding and applications with RGB-D cameras[END_REF]. Furthermore, if depth and colour information can be merged into a single representation, then a complete 3-d representation is possible, in principle. This is clearly desirable, because colour and texture data are essential to many other aspects of scene-understanding, such as identification and tracking [START_REF] Alahi | Cascade of descriptors to detect and track objects across any network of cameras[END_REF].

There are two major obstacles to the construction of a complete scene representation, from a multi-modal camera network. Firstly, typical depth sensors are unable to capture rgb data [START_REF]SR4000 Swiss Ranger manual[END_REF]. This means that the depth and colour cameras will have different viewpoints, and so the raw data are inconsistent. Secondly, typical tof and rgb cameras have limited fields of view, and so the depth and colour data are incomplete. This paper addresses both of these problems, by showing how to estimate the geometric relationships in a multi-view, multi-modal camera network. This task will be called cross-calibration.

In order to constrain the problem, two practical constraints are imposed from the outset. Firstly, the system will be based on time-of-flight (tof) cameras, in conjunction with ordinary rgb cameras. The tof cameras are compact, can be properly synchronized, and are industrially specified, e.g., [START_REF]SR4000 Swiss Ranger manual[END_REF]. Secondly, a modular network of tof+rgb units is required. This is so that individual units can be added or removed, in order to optimize the scene-coverage.

Overview

Time-of-flight cameras can, in principle, be geometrically calibrated by standard methods [START_REF] Hansard | Automatic detection of calibration grids in time-of-flight images[END_REF]. This means that each pixel records an estimate of the scene-distance (range) along the corresponding ray. The 3-d structure of a scene can also be reconstructed from two or more ordinary images, via the parallax (e.g. binocular disparity) between corresponding image points. There are many advantages to be gained by combining the range and parallax data. Most obviously, each point in a parallax-based reconstruction can be mapped back into the original images, from which colour and texture can be obtained. Parallax-based reconstructions are, however, difficult to obtain, owing to the difficulty of putting the image points into correspondence. Indeed, it may be impossible to find any correspondences in untextured regions. Furthermore, if a Euclidean reconstruction is required, then the cameras must be calibrated. The accuracy of the resulting reconstruction will also tend to decrease with the distance of the scene from the cameras [START_REF] Verri | Absolute depth estimate in stereopsis[END_REF].

The range data, on the other hand, are often corrupted by noise and surface-scattering. The spatial resolution of current tof sensors is relatively low, the depth-range is limited, and the luminance signal may be unusable for rendering and for classical image processing. It should also be recalled that tof cameras, of the type used here, cannot be used in outdoor lighting conditions. These considerations lead to the idea of a mixed colour and time-of-flight system, as described in [START_REF] Lindner | Time-of-flight sensor calibration for accurate range sensing[END_REF]. Such a system could, in principle, be used to make high-resolution Euclidean reconstructions, including photometric information [START_REF] Kolb | Time-of-flight cameras in computer graphics[END_REF][START_REF] Gandhi | High-resolution depth maps based on TOF-stereo fusion[END_REF].

In order to make full use of a mixed range/parallax system, it is necessary to find the exact geometric relationship between the different devices. In particular, the reprojection of the tof data, into the colour images, must be obtained. This paper is concerned with the estimation of these geometric relationships. Specifically, the aim is to align the range and parallax reconstructions, by a suitable 3-d transformation.

Previous work

Multi-view depth and colour camera-networks, of the kind used here, produce data-streams that are subject to a variety of geometric relationships [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. These relationships depend on the calibration state, relative orientation, and fields of view of the different cameras. It follows that a variety of calibration strategies can be adopted. These are discussed below, with reference to the literature, and contrasted with the approach presented here.

Perhaps the simplest way to combine rgb and tof data is to perform an essentially 2-d registration between the images and depth maps, as reviewed in [START_REF] Nair | A survey on time-of-flight stereo fusion[END_REF]; see also [START_REF] Huhle | Integrating 3D time-of-flight camera data and high resolution images for 3DTV applications[END_REF][START_REF] Bartczak | Dense depth maps from low resolution time-offlight depth and high resolution color views[END_REF][START_REF] Bleiweiss | Fusing time-of-flight depth and color for realtime segmentation and tracking[END_REF]. This 2-d approach, however, can only provide an instantaneous solution, because changes in the scene-structure produce corresponding changes in the image-to-image mapping. Moreover, owing to the different viewpoints, a complete registration will usually be impossible. If the depth camera also produces a reliable intensity image, then photoconsistency can be used as a 3-d calibration criterion. For example, Beder et al. [START_REF] Beder | Photoconsistent relative pose estimation between a PMD 2D3D-camera and multiple intensity cameras[END_REF] (see also [START_REF] Schiller | Calibration of a PMD camera using a planar calibration object together with a multi-camera setup[END_REF][START_REF] Koch | MixIn3D: 3D mixed reality with ToF-camera[END_REF]) reproject the intensities of the depth data into the colour images, and optimize the camera parameters with respect to the photo-consistency.

Zhu et al. [START_REF] Zhu | Reliability fusion of timeof-flight depth and stereo geometry for high quality depth maps[END_REF] (see also [START_REF] Zhu | Fusion of time-of-flight depth and stereo for high accuracy depth maps[END_REF]) present a sensor-fusion framework for the integration of tof depth and binocular disparity information. This method assumes that a dense disparity map is being computed on-line, which is not required by the method presented in this paper. Furthermore, the geometric calibration method [START_REF] Zhu | Reliability fusion of timeof-flight depth and stereo geometry for high quality depth maps[END_REF] requires manual identification of corresponding points, and is based on a weak perspective camera model. In contrast, our method is automatic, and is based on the more appropriate perspective camera model. However, [START_REF] Zhu | Reliability fusion of timeof-flight depth and stereo geometry for high quality depth maps[END_REF] is complementary to our method, in the sense that their sensor-fusion framework (along with a dense stereo-matcher) could be combined with the projective calibration method described below. Wang and Jia [START_REF] Wang | A fusion framework of stereo vision and kinect for high-quality dense depth maps[END_REF] describe a related sensor-fusion framework for Kinect (rather than tof) depth data and colour images.

Another approach to the multi-modal calibration problem is to apply standard methods, as far as possible, to the depth cameras. Wu et al. [START_REF] Wu | Improved 3D depth image estimation algorithm for visual camera[END_REF] describe an example of this approach. Lindner et al. [START_REF] Lindner | Time-of-flight sensor calibration for accurate range sensing[END_REF] analyze the applicability of standard methods to tof cameras, as well as characterizing the accuracy of the depth data. Mure-Dubois and Hügli [START_REF] Mure-Dubois | Fusion of time-of-flight camera point clouds[END_REF] describe the Euclidean alignment of multiple tof point clouds, having calibrated the cameras by standard methods.

Silva et al. [START_REF] Silva | Camera Calibration using a Color-Depth Camera: Points and Lines Based DLT including Radial Distortion[END_REF] describe a cross-calibration methodology that is based on the identification of 3-d lines in the tof data, which are then projected to corresponding 2-d lines in the rgb images. This method does not require a chequerboard or other calibration pattern; it does, however, require the existence and detection of straight depth-edges throughout the scene. The approach of Silva et al. also involves a non-trivial correspondence problem, which in turn influences the calibration accuracy. Our method uses a standard chequerboard pattern, with a known number of vertices, for which the correspondence problem is relatively straightforward.

Zhang and Zhang [START_REF] Zhang | Calibration between depth and color sensors for commodity depth cameras[END_REF] present cross-calibration methodology that is based on plane constraints, as given in [START_REF] Zhang | Flexible camera calibration by viewing a plane from unknown orientations[END_REF]. This has the advantage of not requiring 2-d features to be detected in the (low resolution) tof images. However, this method cannot address the crucial issue of lens distortion, which is considerable in typical tof cameras [START_REF] Kim | Design and calibration of a multi-view TOF sensor fusion system[END_REF]. A related Kinect-based calibration system is described by Herrera et al. [START_REF] Herrera | Joint Depth and Color Camera Calibration with Distortion Correction[END_REF], again using the plane-based method of [START_REF] Zhang | Flexible camera calibration by viewing a plane from unknown orientations[END_REF]. Herrera et al. give a careful analysis of the Kinect intrinsic parameters, including lens and depth distortion. The latter is analyzed in more detail by Teichman et al. [START_REF] Teichman | Unsupervised intrinsic calibration of depth sensors via slam[END_REF]. Our method does require features to be detected in the tof images, but this also makes it straightforward to estimate the lens parameters, using standard techniques.

Mikhelson et al. [START_REF] Mikhelsona | Automatic, fast, online calibration between depth and color cameras[END_REF] describe an automatic method for registering a Euclidean point-cloud (obtained from a Kinect device) to its 2-d imageprojections. This method, like ours, is based on a chequerboard target. However, Mikhelson et al. perform Euclidean 2-d/3-d registration, in contrast to the more general projective 3-d/3-d registration that is described below. Finally, there are methods that perform 3-d/3-d registration of dense data, subject to pointwise adjustments [START_REF] Cui | Algorithms for 3d shape scanning with a depth camera[END_REF]. This strategy can achieve very close registrations, but introduces a more complex optimization problem, which is not fully compatible with the standard calibration pipeline.

Paper organization and contributions

This paper is organized as follows. Section 2.1 briefly reviews some standard material on projective reconstruction, while section 2.2 describes the representation of range data in the present work. The chief contributions of the subsequent sections are as follows: Section 2.3 describes a point-based method that maps a classical multi-view reconstruction (projective or Euclidean) of the scene onto the corresponding tof representation. The data are obtained from a tof+2rgb system, as shown in figure 1. This does not require the colour cameras to be calibrated (although it is necessary to correct for lens distortion). It is established that this model includes a projective-linear approximation of the systematic tof depth-error.

Section 3 addresses the problem of multi-system alignment, which is necessary for complete scene-coverage. It is shown that this can be achieved in a way that is compatible with the individual tof+2rgb calibrations. The complete cross-calibration pipeline, given a collection of chequerboard images, is fully automatic.

Section 4 contains a detailed evaluation of these methods, using several large data-sets, captured by three tof+2rgb systems (i.e. a nine-camera network). In particular, section 4.1 extends the usual concepts of reprojection error [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] to the multi-modal case. Section 4.2 then introduces a new metric for mixed tof/rgb systems, which measures instantaneous sensor noise, as well as calibration error. The appropriateness of the 3-d homography transformation, as opposed to a similarity transformation, is tested in section 4.3. Section 4.4 discusses possible applications of these systems, including some real 3-d reconstruction examples. Conclusions and future directions are discussed in section 5.

The system presented here is based on the approach introduced by Hansard et al. [START_REF] Hansard | Projective Alignment of Range and Parallax Data[END_REF][START_REF] Hansard | Time-of-Flight Cameras: Principles, Methods and Applications[END_REF]. The earlier work has been improved, and extended to the case of multiple tof and colour cameras. In addition, a new evaluation methodology has been developed, as described above. The automatic detection of calibration targets in the low-resolution tof images, which is a pre-requisite for the methods described here, was developed in a separate paper [START_REF] Hansard | Automatic detection of calibration grids in time-of-flight images[END_REF].

Cross-calibration

This section describes the theory of projective alignment, using the following notation. Bold type will be used for vectors and matrices. In particular, points P, Q and planes U, V in the 3-d scene will be represented by column-vectors of homogeneous coordinates, e.g.

P = P △ P 4 and U = U △ U 4 (1)
where

P △ = (P 1 , P 2 , P 3 ) ⊤ and U △ = (U 1 , U 2 , U 3 ) ⊤ .
The homogeneous coordinates are defined up to a non-zero scaling; for example, P ≃ (P △ /P 4 , 1) ⊤ .

In particular, if P 4 = 1, then P △ contains the ordinary space coordinates of the point P. Furthermore, if |U △ | = 1, then U 4 is the signed perpendicular distance of the plane U from the origin, and U △ is the unit normal. The point P is on the plane U if U ⊤ P = 0. The cross product u × v is often expressed as (u) × v, where (u) × is a 3 × 3 antisymmetric matrix. The column-vector of N zeros is written 0 N . Projective cameras are represented by 3 × 4 matrices. For example, the range projection is

q ≃ CQ where C = A 3×3 | b 3×1 (2)
is a block-decomposition of the 3×4 camera matrix. The left and right colour cameras C ℓ and C r are similarly defined, e.g. p ℓ ≃ C ℓ P. Table 1 summarizes the geometric objects that will be aligned.

Observed Reconstructed Points

Points Planes Points and planes in the two systems are related by the unknown 4 × 4 space-homography H, so that

Binocular C ℓ ,C r p ℓ , p r P U Range C (q, ρ) Q V
Q ≃ HP and V ≃ H -⊤ U. (3)
This model encompasses all rigid, similarity and affine transformations in 3-d.

It preserves collinearity and flatness, and is linear in homogeneous coordinates. Note that, in the reprojection process, H can be interpreted as a modification of the camera matrices. For example, p ℓ ≃ C ℓ H -1 Q, where H -1 Q ≃ P is the point that would theoretically be reconstructed by triangulation. The 4 × 4 homographies H include, as special cases, the rigid transformations that would align a fully-calibrated Euclidean stereo-reconstruction to the tof measurements. There are two motivations for the generalization. Firstly, it allows uncalibrated binocular reconstructions to be used, as described above. Secondly, it has been shown elsewhere [START_REF] Kolb | Time-of-flight cameras in computer graphics[END_REF][START_REF] Zhu | Reliability fusion of timeof-flight depth and stereo geometry for high quality depth maps[END_REF] that tof data are subject to systematic depth biases and nonlinear distortions. These are difficult to correct, owing to the lack of a complete parametric model, and to their dependence on the camera settings (e.g. integration time). Nonetheless, the homography model ( 3) effectively includes a projective-linear approximation of the depth distortion, which is fitted along with the other transformation parameters. This model is quite powerful: it includes rational depth-distortions, with varying parameters, across each bundle of rays. For example, the two-parameter inverse disparity calibration, as used with Kinect devices [START_REF] Herrera | Joint Depth and Color Camera Calibration with Distortion Correction[END_REF], is a special case of the homography model described here. Indeed, even if the rgb cameras are fully calibrated, the 3-d homographies are needed to account for depth-distortions and residual reconstruction errors, as demonstrated in [START_REF] Hansard | Projective Alignment of Range and Parallax Data[END_REF]. This issue will be explored in section 4.3, below.

Parallax-based reconstruction

A projective reconstruction of the scene can be obtained from matched points p ℓk and p rk , together with the fundamental matrix F, where p ⊤ rk F p ℓk = 0. The fundamental matrix can be estimated automatically, using the wellestablished ransac method. The camera matrices can then be determined, up to a four-parameter projective ambiguity [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. In particular, from F and the epipole e r , the cameras can be defined as

C ℓ ≃ I | 0 3 ) and C r ≃ (e r ) × F + e r g ⊤ γe r . (4)
where γ = 0 and g = (g 1 , g 2 , g 3 ) ⊤ can be used to bring the cameras into a plausible form. This makes it easier to visualize the projective reconstruction and, more importantly, can improve the numerical conditioning of subsequent procedures.

Range-based reconstruction

The tof camera C provides the range (i.e. radial distance) ρ of each scenepoint from the camera-centre, as well as the associated image-coordinates q = (x, y, 1). The back-projection of this point into the scene is

Q △ = A -1 (ρ/α) q -b where α = A -1 q . ( 5
)
Hence the point (Q △ , 1) ⊤ is at distance ρ from the optical centre -A -1 b, in the direction A -1 q. The scalar α serves to normalize the direction-vector. This is the standard pinhole model, as used in [START_REF] Beder | A comparison of PMD-cameras and stereo-vision for the task of surface reconstruction using patchlets[END_REF].

The range data are noisy and incomplete, owing to illumination and scattering effects. This means that, given a sparse set of features in the intensity image (of the tof device), it is not advisable to use the back-projected point (5) directly. A better approach is to segment the image of the plane in each tof camera (using the the range and/or intensity data). It is then possible to back-project all of the enclosed points, and to robustly fit a plane V j to the enclosed points Q ij , so that V ⊤ j Q ij ≈ 0 if point i lies on plane j. Now, the back-projection Q π of each sparse feature point q can be obtained by intersecting the corresponding ray with the plane V, so that the new range estimate ρ π is

ρ π = V ⊤ △ A -1 b -V 4 (1/α) V ⊤ △ A -1 q (6)
where |V 4 | is the distance of the plane to the camera centre, and V △ is the unit-normal of the range plane. The new point Q π is obtained by substituting ρ π into (5).

The choice of plane-fitting method is affected by two issues. Firstly, there may be very severe outliers in the data, due to the photometric and geometric errors in the depth-estimation process. Secondly, the noise-model should be based on the pinhole geometry, which means that perturbations occur radially along visual directions, which are not (in general) perpendicular to the observed plane [START_REF] Hebert | 3D measurements from imaging laser radars: How good are they?[END_REF][START_REF] Wang | Comparison of local plane fitting methods for range data[END_REF]. Several plane-fitting methods, both iterative [START_REF] Kanazawa | Reliability of Plane Fitting by Range Sensing[END_REF] and non-iterative [START_REF] Pathak | Revisiting uncertainty analysis for optimum planes extracted from 3D range sensor point-clouds[END_REF] have been proposed for the pinhole model. However, for tof data, the chief problem is the large number of outliers. This means that a ransac-based method is the most effective in this context [START_REF] Hansard | Projective Alignment of Range and Parallax Data[END_REF].

Projective Alignment

It is straightforward to show that the transformation H in (3) could be estimated from five binocular points P k , together with the corresponding range points Q k . This would provide 5 × 3 equations, which determine the 4 × 4 entries of H, subject to an overall projective scaling. It is better, however, to use the 'Direct Linear Transformation' method [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], which fits H to all of the data. This method is based on the fact that if

P ′ = HP (7)
is a perfect match for Q, then µQ = λP ′ , and the scalars λ and µ can be eliminated between pairs of the four implied equations [START_REF] Csurka | Finding the collineation between two projective reconstructions[END_REF]. This results in homogeneous equations as

Q 4 P ′ △ -P ′ 4 Q △ Q △ × P ′ △ = 0 6 . (8)
Note that if P ′ and Q are normalized so that P ′ 4 = 1 and Q 4 = 1, then the magnitude of the top half of ( 8) is simply the distance between the points. Following Förstner [START_REF] Förstner | Uncertainty and projective geometry[END_REF], the left-hand side of ( 8) can be expressed as Q ∧ P ′ where

Q ∧ = Q 4 I 3 -Q △ Q △ × 0 3 (9)
is a 6 × 4 matrix, and Q △ × P △ = Q △ × P △ , as usual. The equations ( 8) can now be written in terms of ( 7) and ( 9) as

Q ∧ HP = 0 6 . (10)
This system of equations is linear in the unknown entries of H, the columns of which can be stacked into the 16 × 1 vector h. The Kronecker product identity vec(X Y Z) = (Z ⊤ ⊗ X) vec(Y) can now be applied, to give

P ⊤ ⊗ Q ∧ h = 0 6 where h = vec H . ( 11
)
If M points are observed on each of N planes, then there are k = 1, . . . , M N observed pairs of points, P k from the projective reconstruction and Q k from the range back-projection. The M N corresponding 6×16 matrices P ⊤ k ⊗ (Q k ) ∧ are stacked together, to give the complete system

   P ⊤ 1 ⊗ Q 1 ∧ . . . P ⊤ M N ⊗ Q M N ∧    h = 0 6M N ( 12
)
subject to the constraint |h| = 1, which excludes the trivial solution h = 0 16 . It is straightforward to obtain an estimate of h from the SVD of the the 6M N × 16 matrix on the left of [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. This solution, which minimizes an algebraic error [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], is the singular vector corresponding to the smallest singular value of the matrix. Note that the point coordinates should be transformed, to ensure that ( 12) is numerically well-conditioned [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. In this case the transformation ensures that k P k△ = 0 3 and 1

M N k |P k△ | = √ 3,
where P k4 = 1. The analogous transformation is applied to the range points

Q k .
The above procedure effectively computes a projective basis in 3-d, which requires five points, no four of which may be co-planar. This gives 5 × 3 numbers, which are equivalent to the 15 parameters of the homogeneous 4 × 4 matrix H. In practice, however, we require full detection of 35 vertices per board, and around 10 boards per homography (depending on visibility constraints). Hence the procedure operates far from any degenerate cases, which would involve fewer than 15 degrees of freedom in total.

The DLT method, in practice, gives a good approximation H DLT of the homography (3). This can be used as a starting-point for the iterative minimization of a more appropriate error measure. In particular, consider the reprojection error in the left image,

E ℓ (C ℓ ) = M N k=1 D C ℓ Q k , p ℓk 2 (13)
where D(p, q) = | p △ /p 3q △ /q 3 |. A 12-parameter minimization of the reprojection error [START_REF] Nair | A survey on time-of-flight stereo fusion[END_REF], starting with the linear estimate

C ℓ ← C ℓ H -1
DLT , is then performed by the Levenberg-Marquardt algorithm [START_REF] Press | Numerical Recipes in C[END_REF]. The result will be the camera matrix C ⋆ ℓ that best reprojects the range data into the left image (C ⋆ r is similarly obtained). The solution, provided that the calibration points adequately covered the scene volume, will remain valid for subsequent depth and range data.

Alternatively, it is possible to minimize the joint reprojection error, defined as the sum of left and right contributions,

E H -1 = E ℓ C ℓ H -1 + E r C r H -1 (14)
over the (inverse) homography H -1 . The 16 parameters are again minimized by the Levenberg-Marquardt algorithm, starting from the DLT solution H -1 DLT . The difference between the separate (13) and joint [START_REF] Huhle | Integrating 3D time-of-flight camera data and high resolution images for 3DTV applications[END_REF] minimizations is that the latter preserves the original epipolar geometry, whereas the former does not. Recall that C ℓ C r , H and F are all defined up to scale, and that F satisfies an additional rank-two constraint [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Hence the underlying parameters can be counted as (12 -1) + (12 -1) = 22 in the separate minimizations, and as (16 -1) = 15 in the joint minimization. The fixed epipolar geometry accounts for the (9 -2) missing parameters in the joint minimization. If F is known to be very accurate (or must be preserved) then the joint minimization ( 14) should be performed. This will also preserve the original binocular triangulation, provided that a projective-invariant method was used [START_REF] Hartley | [END_REF]. However, if minimal reprojection error is the objective, then the cameras should be treated separately. This will lead to a new fundamental matrix

F ⋆ = (e ⋆ r ) × C ⋆ r (C ⋆ ℓ ) +
, where (•) + is the generalized inverse, and C ⋆ ℓ , C ⋆ r are the optimized camera-matrices. The epipole in the right-hand image is obtained from e ⋆ r = C ⋆ r d ⋆ ℓ , where the vector d ⋆ ℓ represents the nullspace

C ⋆ ℓ d ⋆ ℓ = 0 3 .

Multi-System Alignment

The methods described in section 2.3 can be used to calibrate a single tof+2rgb system; the joint calibration of several such systems will now be explained. In this section the notation P i will be used for the binocular coordinates (with respect to the left camera) of a point in the i-th system, and likewise Q i for the tof coordinates of a point in the same system. Hence the i-th tof, left and right rgb cameras (sharing the same physical mounting) have the form

C i ≃ A i | 0 3 ) C ℓi ≃ A ℓi | 0 3 ) and C ri ≃ A ri | b ri ) (15)
where A i and A ℓi contain only intrinsic parameters, whereas A ri also encodes the relative orientation of C ri with respect to C ℓi . Each system has a transformation H -1 i that maps tof points Q i into the corresponding rgb coordinate system of C ℓi . Furthermore, let the 4 × 4 matrix G ij be the transformation from system j, mapping back to system i (between different physical mountings). This matrix, in the binocularly-calibrated case, is a rigid 3-d transformation. However, by analogy with the tof-to-rgb matrices H, each G ij could be a projective transformation in the uncalibrated case; this would allow Euclidean structure to propagate from the tof measurements, across the entire camera network.

The left and right cameras, in all cases, that project a scene-point P j in coordinate system j to image-points p ℓi and p ri in system i are

C ℓij = C ℓi G ij and C rij = C ri G ij . (16)
Note that if a single global coordinate system is chosen to coincide with the k-th rgb system, then a point P k projects via C ℓik and C rik . These two cameras are respectively equal to C ℓi and C ri in [START_REF] Bartczak | Dense depth maps from low resolution time-offlight depth and high resolution color views[END_REF] only when i = k, such that G ij = I in [START_REF] Bleiweiss | Fusing time-of-flight depth and color for realtime segmentation and tracking[END_REF]. A typical three-system configuration is shown in fig. 2. Each ellipse represents a separate system, with system 2 chosen as the reference. The arrows (with camera-labels) show some possible tof-to-rgb projections. For example, a point P 2 ≃ H -1 2 Q 2 in the centre projects directly to rgb view ℓ2 via C ℓ2 , whereas the same point projects to ℓ3 via

C ℓ32 = C ℓ3 G 32 .
The transformation G ij can only be estimated directly if there is a region of common visibility between systems i and j. If this is not the case (as when the systems face each other, such that the front of the calibration board is not simultaneously visible), then G ij can be computed indirectly. For example, G 02 = G 01 G 12 where P 2 = G -1 12 G -1 01 P 0 . Note that, in all cases, the stereo-reconstructed points P are used to estimate these transformations This is because they are always more reliable than the tof points Q, as demonstrated below.

Evaluation

The following sections will describe the accuracy of a nine-camera setup, calibrated by the methods described above. Section 4.1 will evaluate calibration error, whereas section 4.2 will evaluate total error. The former is essentially a fixed function of the estimated camera matrices, for a given scene. The latter also includes the range-noise from the tof cameras, which varies from moment to moment (due to intrinsic noise, changing illumination, and object motion). The importance of this distinction will be discussed. In section 4.3 we analyze the avantages of using homographies to align the data, as opposed to similarity transformations. Finally, in section 4.4 we demonstrate some applications of the complete system.

The setup consists of three rail-mounted tof+2rgb systems, i = 1 . . . 3, as in fig. 2. The stereo baselines are 17cm on average, and the tof cameras are separated by 107cm on average. The rgb images are 1624 × 1224, whereas the Mesa Imaging sr4000 tof images are 176 × 144, with a depth range of 500cm [START_REF]SR4000 Swiss Ranger manual[END_REF]. The three stereo systems are first calibrated by standard methods [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], returning a full Euclidean decomposition of C ℓi and C ri , as well as the associated lens parameters. The lenses of the tof cameras are also calibrated by standard methods [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. This removes some of the radial depth deformation that has been observed in the tof data [START_REF] Cui | Algorithms for 3d shape scanning with a depth camera[END_REF]. The matrices G ij are rigid-body transformations, which are estimated by the SVD method of Arun et al. [START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF]. Projective alignment is preferred for the tof/rgb alignment, for reasons discussed in sections 2 and 10. Hence the transformations H -1 j are 4 × 4 homographies, estimated by the method of section 2.3. Specifically, the DLT solutions were refined by Levenberg-Marquardt minimization of the joint geometric error, as in ( 14).

Calibration Error

The calibration error is measured by first taking tof points Q π j corresponding to chequerboard vertices on the reconstructed calibration plane π j in system j, as described in section 2.2. These can then be projected into a pair of rgb images in system i, so that the geometric image-error

E cal ij = 1 2 E cal ℓij + E cal rij can be computed
, where

E cal ℓij = 1 |π| Q π j D C ℓij H -1 j Q π j , p ℓi (17)
and E cal rij is similarly defined. The function D(•, •) computes the imagedistance between two inhomogenized points, as in [START_REF] Nair | A survey on time-of-flight stereo fusion[END_REF], and the denominator corresponds to the number of vertices on the board, with |π| = 35 in the present experiments. The measure ( 17) can of course be averaged over all images in which the board is visible. The rgb cameras were calibrated, by standard methods, as described above. Subpixel accuracy was obtained, which confirms the accuracy of the camera matrices C ℓi and C ri , as well as the inter-system matrices G ij , where i, j = 1, 2, 3.

For the purpose of evaluation, a new set of tof-vertices Q π i were reconstructed, fitted and reprojected within each system i. This evaluation effectively tests the quality of the H -1 i matrices, by comparing C ℓi H -1 i Q π i to p ℓi , and analogously to points p ri in the other image. Note that all camera and transformations parameters are now fixed; no optimization was performed with respect to the evaluation data. The whole experiment was performed on three large data-sets, from different capture-sessions, and with different camera configurations, labelled A, B and C. Such a configuration leads to one triplet per system as shown in fig. 3. While an image can be frontoparallel for one frame, it may appear very slanted to the other frames. The example of fig. 3 shows a calibration image that is almost fronto-parallel to the tof frame of the central tof+2rgb unit. Table 2 shows that the average reprojection error remains subpixel in all three data-sets. The corresponding error-distributions are shown as histograms in fig. 4. It is also interesting to consider how the above tof-vertices reproject into different systems. The error 3-d error distribution, for a given pixel (either tof or stereo) is highly anisotropic; the direction of the corresponding ray is much more reliable than the distance along it [START_REF] Cui | Algorithms for 3d shape scanning with a depth camera[END_REF]. In practice, the different tof+2rgb systems are distributed around the edge of the room, looking inwards. Hence a given system is likely to be seen 'from the side' by at least one other system. This means that any large depth errors in the first system will not cancel-out in the re-projection to the other systems. This effect is seen clearly in figure 5, which shows calibration errors of up to several pixels, from one system to the others.

Set Mean Median

Total Error

The calibration error, as reported in the preceding section, is the natural way to evaluate the estimated cameras and homographies. It is not, however, truly representative of the 'live' performance of the complete setup. This is because the calibration error uses each estimated plane π j to replace all vertices Q j with the fitted versions Q π j . In general, however, no surface model is available, and so the raw points Q j must be used as input to segmentation, meshing and rendering processes.

The total error, which combines the calibration and range errors, can be measured as follows. The i-th rgb views of plane π j are related to the tof image-points q j by the 2-d transfer homographies T ℓij and T rij , where These 3×3 matrices can be estimated to subpixel accuracy, by using the DLT algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] to obtain initial homographies, and then applying area-based alignment [START_REF] Evangelidis | Parametric image alignment using enhanced correlation coefficient maximisation[END_REF] to produce the final estimates T ℓij and T rij .

p ℓi ≃ T ℓij q j and p ri ≃ T rij q j . (18)
Let Π j be the 2-d hull (i.e. bounding-polygon) of plane π j as it appears in the tof image. Any pixel q j in the hull (including the original calibration vertices) can now be re-projected to the i-th rgb views via the 3-d point Q j , or transferred directly by T ℓij and T rij in [START_REF] Schiller | Calibration of a PMD camera using a planar calibration object together with a multi-camera setup[END_REF], as shown in figure 6. In other words, we are able to isolate the image-to-image error that is incurred by mapping via the 3-d reconstruction, in relation to the direct 2-d to 2-d mapping defined by T ℓij and T rij .

The total error is the average difference between the reprojections and the transfers,

E tot ij = 1 2 E tot ℓij + E tot rij
, where

E tot ℓij = 1 |Π j | q j ∈Π j D C ℓij H -1 j Q j , T ℓij q j (19)

and E tot

rij is similarly defined. The view-dependent denominator |Π j | ≫ |π| is the number of pixels in the hull Π j . Hence E tot ij is the total error, including both calibration errors and range-noise, of tof plane π j as it appears in the i-th rgb cameras.

If the rgb cameras are not too far from the tof camera, then the range errors tend to be cancelled in the reprojection. This is evident from table 3, in which the total errors, on average, are not much greater than the calibration errors in table 2. It is, however, clear from fig. 7, that the tails of the total error distributions are greatly increased by outliers in the tof data-stream. Although errors up to 3 pixels are shown in distributions for the sake of clarity, maximum errors of more than five pixels are common (see table 3). [START_REF] Koch | MixIn3D: 3D mixed reality with ToF-camera[END_REF], measured by projecting the raw tof points Q i to the left and right rgb images (1624 × 1224) of the respective systems i = 1, 2, 3. These figures characterize the raw data that is produced by the live system. This data includes all outliers in the tof data-stream, which gives rise to at least one error of 20 pixels here.

Set Mean Median

When the raw tof points are reprojected to a different system, a high total error is expected, because the sensor's range-noise may be viewed 'from the side'. Indeed, fig. 8 shows that a substantial proportion of the tof points reproject with total errors in excess of ten pixels. These kind of gross errors are characteristic of the tof data, owing to the inevitable presence of absorbing and scattering surfaces in a typical scene. In fact, calibration error and total error are both influenced by surface orientation. For example, the higher errors in Set B of the data (see tables 2 and 3) are due to the presence of some very slanted boards, in both the fitting and evaluation data-sets; this causes two problems, as follows. Firstly, the images of these boards are very foreshortened in some systems. This leads to increased calibration error, because the vertices are hard to detect in the corresponding tof images. Secondly, the strength of the reflected IR signal is reduced whenever the surface is oblique to a given tof camera. If the surface is also absorbent (like the black squares of the board), then the total error is greatly increased, due to the combined effects of scattering and absorption.

It is possible to understand these results more fully by examining the distribution of the total error across individual boards. Figure 9 shows the distribution for a board reprojected to same/different systems (i.e. part of the data from figs. 7 & 8). There is a relatively smooth gradient of error across the board, which is attributable to errors in the fitting of plane π j , and in the estimation of the camera parameters. In addition to these effects, it is clear that the gross errors are correlated with the black squares of the board, which reflect too little of the tof signal. For instance, in the bottom example of figure 9, the mean reprojection error associated with the black squares is 10.05 pixels, whereas the white squares are associated with a mean error of 3.55 pixels. The effect is particularly noticeable, as expected from the histograms, when reprojecting to different camera systems.

Comparison of homography and similarity transformations

It was argued in section 2 that the 3-d homography model is an appropriate way to compensate for miscalibrations of the tof and stereo systems (and indeed it allows the stereo system to be left uncalibrated, if preferred [START_REF] Hansard | Projective Alignment of Range and Parallax Data[END_REF]). This claim is tested in the following experiments, using data images from two new data sets.

The most obvious alternative to the homography is a rigid motion and scaling; i.e. a 3-d similarity transformation. Specifically, the 4 × 4 matrix H in (3) can be replaced by the similarity transformation

S = σ R t 0 1 ( 20
)
where σ is a (positive) scalar, R is a 3 × 3 rotation matrix, and t is a 3 × 1 translation vector. Hence there are only seven degrees of freedom, in contrast to the fifteen of the homography. Geometrically, the similarity model can be interpreted as a rigid transformation between the tof and stereo systems, where the baseline distance of the latter is unknown.

In principle, a homography should always result in equal or lower error, because it includes the similarity transformation as a special case. In particular, a homography can always be written as H = SD where the deformation D is composed of an affinity and an elation [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. However, there are two practical issues to consider. Firstly, the quality of the actual estimate in each case; and secondly, the possible danger of over-fitting with the homography.

It was argued in section 2 that the 3-d homography model is an appropriate way to compensate for miscalibrations of the tof and stereo systems (and also allows the stereo system to be left uncalibrated, if preferred [START_REF] Hansard | Projective Alignment of Range and Parallax Data[END_REF]). In particular the deformation D can help to account for locally linear depth bias in the tof camera [START_REF] Foix | Lock-in Time-of-Flight (ToF) Cameras: A Survey[END_REF]. Hence the additional eight degrees of freedom in the homography model allow for an overall approximation to the per-pixel corrections performed by Cui et al. [START_REF] Cui | Algorithms for 3d shape scanning with a depth camera[END_REF].

The procedure to estimate the optimal similarity (20) is analogous to that used to estimate the homography in section 2. Instead of DLT, an initial estimate is obtained from the Procrustes algorithm [START_REF] Schönemann | Fitting one matrix to another under choice of a central dilation and a rigid motion[END_REF][START_REF] Horn | Closed-form solution of absolute orientation using orthonormal matrices[END_REF], using three or more correspondences. The reprojection error is then minimized using the Levenberg-Marquardt procedure, as with the homography. Note that the rotation R in ( 20) is appropriately parameterized by the Rodrigues formula.

The results of these experiments are shown in figure 10. In every tof+2rgb system, and in both capture sessions, the homography results in lower error than the similarity. The mean reprojections errors (in pixels) were 0.22 vs. 0.65, and 0.46 vs. 1.06, for the respective capture sessions. The generally higher figures from the second set are due to a broader distribution of board-poses in the evaluation data-set.

These experiments show that the 3-d homography transformation is a suitable model for tof/rgb alignment, as argued in section 2. It may also be noted that the homography is effectively easier to estimate than the similarity, as there is no need to maintain the orthogonality of R during the final minimization.

Cross-calibration for 3-d reconstruction and rendering

The proposed cross-calibration methodology also allows a fully-textured 3-d model to be segmented and rendered. This is because the depth-data are automatically assigned to pixels in the rgb images. Figure 11 shows reconstruction instances with and without texture from a 360 • reconstruction, obtained from the cross-calibration of four tof+2rgb systems. Meshing was performed by the standard Poisson reconstruction method [START_REF] Kazhdan | Poisson surface reconstruction[END_REF], followed by simple triangle-based texture mapping.

More dense reconstruction can be achieved by exploiting the full resolution of stereo cameras, as described elsewhere [START_REF] Gandhi | High-resolution depth maps based on TOF-stereo fusion[END_REF]. Figure 12 shows the difference between the raw reconstruction obtained from the cross-calibration of a tof+2rgb system, and that obtained after using the tof data in a subsequent stereo algorithm [START_REF] Gandhi | High-resolution depth maps based on TOF-stereo fusion[END_REF]. A high-resolution depth map is produced, making full use of the high-resolution rgb cameras.

Conclusion

It has been shown that 3-d projective transformations can be used to cross-calibrate a system of tof and stereo reconstructions. A practical method for computing these transformations, based on geometric principles, has been introduced. This calibration procedure has been extended to a The raw mesh, with textured edges, obtained directly from our tof+2rgb cross-calibration method. Right: The textured high-density point cloud, obtained after applying the tof-stereo fusion method of [START_REF] Gandhi | High-resolution depth maps based on TOF-stereo fusion[END_REF] to the cross-calibration results.

nine-camera network of six rgb and three tof cameras, and evaluated in detail.

A clear distinction has been made between calibration error, which is due to imperfect camera and image models, versus total error, which incorporates tof-specific noise and biases. It has been shown that these errors can be separated geometrically, and used to characterize the performance of a tof/rgb camera network. The overall performance of the system has been visualized in segmented 360 • 3-d reconstructions. This type of reconstruction is challenging to compute, and shows that the system is more than adequate as a basis for scene-segmentation tasks.

The accuracy of the method presented here is somewhat limited by the relatively poor localization of points, even when spatially interpolated, in the tof images. The detection of standard chequerboard patterns, in tof images, has been discussed elsewhere [START_REF] Hansard | Automatic detection of calibration grids in time-of-flight images[END_REF]. The design of more convenient calibration patterns, for tof cameras, would be a worthwhile direction for future research. Meanwhile, in mitigation, the spatial resolution of tof cameras will continue to increase.

Future work should also consider the distribution of range-errors in 3-d, and how this can be used to design custom meshing and surface reconstruction algorithms for time-of-flight data. This, in conjunction with the rgb textures provided by the present method, would lead to a comprehensive approach to 3-d reconstruction, rendering and scene-understanding.

Figure 1 :

 1 Figure 1: A single tof+2rgb system (left), comprising a time-of-flight camera in the centre, plus a pair of ordinary colour cameras. This paper addresses the problem of simultaneously cross-calibrating many such systems (right), as a foundation for scene understanding

Figure 2 :

 2 Figure 2: Example of a three tof+2rgb setup, with tof cameras labelled 1,2,3; cf. fig. 1.Each ellipse represents a separate system, with system 2 chosen as the reference. The arrows (with camera-labels) show some possible tof-to-rgb projections. For example, a point P 2 ≃ H -1 2 Q 2 in the centre projects directly to rgb view ℓ2 via C ℓ2 , whereas the same point projects to ℓ3 via C ℓ32 = C ℓ3 G 32 .

  17), measured by projecting the fitted tof vertices Q π i to the left and right rgb images (1624 × 1224) of the respective systems i = 1, 2, 3. The experiment was repeated three times (A-C), with different camera configurations. Each statistic was computed from the left and right rgb-reprojections of 35 vertices in 7 views of the board (total number of 2-d points, per data-set: 3 × 2 × 35 × 7 = 1470).

Figure 3 :

 3 Figure 3: Calibration images from synchronized captures. Each row corresponds to the image triplet of a tof+2rgb system, while the range images are shown as well. Grayscale images of tof sensors have undergone lens undistortion and contrast enhancement.

Figure 4 :

 4 Figure 4: Calibration error (17), measured by projecting the fitted tof points Q π i to the left and right rgb images (1624 × 1224) of the respective systems i = 1, 2, 3. The experiment was repeated three times (A-C), with different camera configurations. Each histogram contains 1470 points, as in table 2.

Figure 5 :

 5 Figure 5: Calibration error (17), measured by projecting the fitted tof points Q π 2 to the left and right rgb images (1624 × 1224) of two different systems, i = 1, 3. Each histogram combines left-camera and right-camera measurements from 7 views of the calibration board.

Figure 6 :

 6 Figure 6: The 2-d transfer homography T ℓ describes the chessboard mapping from the tof frame to the left frame. In principle, this mapping should be equivalent to a mapping via the range data, i.e. from the tof frame to the measured 3-d plane, and back to the left frame.

Figure 7 :

 7 Figure 7: Total error (19), measured by projecting the raw tof points Q i to the left and right rgb images (1624 × 1224) of the respective systems i = 1, 2, 3. These distributions have longer and heavier tails than those of the corresponding calibration errors, shown in fig. 4.

Figure 8 :

 8 Figure 8: Total error[START_REF] Koch | MixIn3D: 3D mixed reality with ToF-camera[END_REF], measured by projecting the fitted tof points Q π 2 to the left and right rgb images (1624 × 1224) of two different systems, i = 1, 3. The range errors are emphasized by the difference in viewpoints between the two systems. Average error is now around five pixels, and the noisiest tof points reproject with tens of pixels of error.

Figure 9 :

 9 Figure9: 3-d tof pixels on an example calibration board, projected into an RGB camera in the same (top) and in a different (bottom) system. Each pixel is colour-coded according to the relative total-error[START_REF] Koch | MixIn3D: 3D mixed reality with ToF-camera[END_REF] in each case (larger errors are shown in yellow). Black crosses mark the detected vertices in the rgb image. The black squares of the board are associated with larger reprojection errors, particularly after reprojecting to a different system (mean error in the bottom example is 10.05px for black squares, and 3.55px for white).

Figure 10 :

 10 Figure 10: Calibration error (17), measured by projecting the fitted tof points Q π i to the left and right rgb images (1624 × 1224) of the respective systems i = 1, 2, 3. The tof/stereo alignment was computed via a 3-d homography (top), or a similarity transformation (bottom). It is clear that the former is superior, as described in section 4.3. The left and right 2 × 3 blocks represent different capture-sessions, as described in the text.

Figure 11 :

 11 Figure 11: Examples of a segmented 360 • reconstruction, after Poisson meshing, with and without rgb texture.

Figure 12 :

 12 Figure 12: Cross-calibration for 3-d reconstruction. Left:The raw mesh, with textured edges, obtained directly from our tof+2rgb cross-calibration method. Right: The textured high-density point cloud, obtained after applying the tof-stereo fusion method of[START_REF] Gandhi | High-resolution depth maps based on TOF-stereo fusion[END_REF] to the cross-calibration results.

Table 1 :

 1 Summary of notations in the left, right and range systems.

Table 2 :

 2 Calibration error (

				Max Count
	A	0.59	0.52	1.82 1470
	B	0.72	0.58	4.86 1470
	C	0.45	0.40	1.48 1470

Table 3 :

 3 Total error

2 = 6 interdependent constraints per point. It is convenient to write these

This work has received funding from Agence Nationale de la Recherche under the MIXCAM project number ANR-13-BS02-0010-01.


Mario Sanz Lopez

Jérémie Dequidt

Erwan Kerrien

Christian Duriez

Marie-Odile Berger

Stéphane Cotin

Mario Sanz-Lopez

Stephane Cotin

Testbed for assessing the accuracy of interventional radiology simulations

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Interventional radiology procedures rely on the navigation of a catheter inside the blood vessels through fluoroscopic image guidance. When the tip of the catheter is close to the location of the disease, dedicated tools (coils, stents, balloon. . . ) can be deployed or drugs can be delivered in order to diagnose or treat the patient. To control the motion of a catheter within the vascular network, the radiologist can only push, pull or twist the proximal end of the device. She/he also can insert a stiffer device, named guidewire, that will help reaching specific locations. Since such devices are constrained inside the patient's vasculature, the motion of the devices results from the displacement inputs of the radiologist as well as contact forces between the devices and the arterial walls. These contact forces are particularly important since the motion of the devices is greatly affected by dynamic effects such as friction or sliding which can significantly change the trajectory followed by the devices. Even if the friction between the arterial wall and the tools may be negligible in certain circumstances (for instance slow motion, simple path to reach the targeted area), its effects are significant in tortuous vasculature or aneurysms. The figure 1 illustrates that friction implies high frequency motion for the deployment of a coil (thin coated platinum thread) in an aneurysm. Therefore, accurate models for the tools and for the interactions between tools and the blood vessel are mandatory. However, the quantitative evaluation of the accuracy is hard to perform due to the complexity of the environment: while imaging techniques allow to capture the geometry of the anatomical structures, the bio-mechanical parameters are hard to measure / estimate. Some previous work has proposed strategies to estimate the accuracy of their simulations: Dequidt et. al. [START_REF] Dequidt | Interactive simulation of embolization coils: Modeling and experimental validation[END_REF] use 3D rotational angiography to evaluate the simulation of a coil deployment in a collision-free environment. Luboz et. al. [START_REF] Luboz | Simulation of endovascular guidewire behaviour and experimental validation[END_REF] compare the position of simulated guidewires and actual guidewires inserted in a vascular phantom. The comparison is however performed on equilibrium positions and on a almost linear trajectory which reduces the influence of friction. Finally, Alderliesten et. al [START_REF] Alderliesten | Modeling friction, intrinsic curvature, and rotation of guide wires for simulation of minimally invasive vascular interventions[END_REF] have conducted a thorough quantitative experimental validation of a quasi-static model of guidewires: 3D comparison through 3DRA is performed with several deployment scenarios. However the main limitation of the approach is that only static friction is considered and the friction parameter is empirically set. Even if the friction between the arterial wall and the tools may be negligible in certain circumstances, its effects are significant in tortuous vasculature or aneurysms even if special care is taken by the radiologist with slow and delicate manipulations. For instance, friction implies high frequency motion for the deployment of a coil (thin coated platinum thread) in an aneurysm. To capture this effect on navigation scenarios, our objective is to design a testbed with software support that is able to record dynamic and non-smooth phenomena (i.e high speed changes). This paper describes an experimental testbed suited for navigation scenarios that include motion tracking of the tool but also image tracking through high-speed stereoscopic cameras.

Validation testbed

The testbed proposed in this paper is designed to provide ground-truth data to compare against simulation. To do so, the testbed is composed of two components: (1) a tracking device is built to accurately track the input motion of the surgeon navigating the radiology tools, (2) high-speed stereoscopic cameras coupled with segmentation and reconstruction provide the position of the tool in a vascular silicon phantom.

Motion Tracking:

We have designed a new motion tracking interface, based on an optical sensor from a mouse that detects motion with a 1200 dpi resolution, a standard HID USB interface and a radio wireless connexion to transmit the raw data without cluttering the operating field. A new package has been built Fig. 1: Recorded sequence of a coil deployed at low and constant speed where image acquisition is performed at 15 fps through 3D Rotational Angiography. The sequence shows the effect of friction and then sliding as the coil completely changes its configuration between frames 3 and 4 (67 msecs between each frame).

for the sensor allowing an easy handle, standard surgical attach systems on both sides and removable tunnels adapted to each size of catheters and coils. The main issues with the use of a mouse sensor relies in sensing the cylinder shape of the catheter instead of a flat surface, the placement of the circuit in a different embodiment, as well as the existence of a wireless, potentially noisy link rise questions about its reliability. Moreover, the data provided is differential, which could lead to a cumulative input error. Thus, a precision test has been conducted considering both local (resolution) and absolute accuracy using a calibrated test bench.

This test bench is composed of a motorized spool that pulls over one meter of bare electric wire, with a negligible elongation factor and a constant diameter of 1.54mm (close to interventional radiology tools). While the electric wire diameter is greater that the actual interventional radiology tools, it was chosen due to its non-elongation property and we use dedicated diameter-dependent removable tunnels that allow to place any thin wire in the focal plane of the sensor in order to maximize the motion tracking. Moreover the wire is not insulated as it does not provide a significant difference for the tracking. The motion of the spool is tracked using the encoder of a ball mouse and controlled through a microcontroller connected to power interfaces driving two motors at both sides of the cylinder; the grip between the encoder wheel and the spool is ensured by choosing a silicone rubber solid cylinder for the spool, with typical friction coefficient larger than 1. Under these conditions, the mouse sent exactly 1260 counts, which gave a final system absolute error of ±396µm, that can be considered negligible. A landmark was set on the wire in order to compare after each test that the mark was perfectly aligned with the goal position of 1 meter of wire pulled.

The laser mouse used as sensor makes data easy to receive, but it needs to be fetched before the Operating System applies transfer functions and denaturalizes the raw recovered data, this was done with the libpointing library [START_REF] Casiez | No more bricolage!: Methods and tools to characterize, replicate and compare pointing transfer functions[END_REF]. 6 tests were conducted at various speeds and the distance recorded by the sensor was stored. For 1 meter of wire pulled with a speed of (0.0487, 0.0967, 0.145, 0.192, 0.242 and 0.290 m/s) the sensed distance was (1051.200, 1006.284, 1009.480, 1064.662, 1060.026 and 991.4255 mm) which leads to an average error of 3.3%.

High-speed Stereoscopic Tracking:

The system was centered on two JAI-Pulnix TM-6740CL cameras configured to acquire 640x480 grayscale images at 198 frames per second. Synchronization was performed thanks to a trigger (CC320 Machine Vision Trigger Timing Controller, Gardasoft). Two Super Cool-Lite 9 (interFit Lighting) projectors were used to reduce the flickering observed with neon lights. Both cameras were placed to respectively look down and sideways (see Fig. 2, left). The stereo pair of cameras were calibrated using chessboard detection and calibration by OpenCV.

The vasculature was given by a phantom made of silicon (right internal carotid artery rigid model with three aneurysms, by Elastrat). Silicon is transparent (see Fig. 2, middle) but induces refraction that significantly perturbs the epipolar geometry of the stereoscopic cameras and forbids using traditional 3D reconstruction based on the triangulation of visual rays [START_REF] Chari | Multiple-view geometry of the refractive plane[END_REF]. The plane equations of both refractive faces were given by computing the associated homography thanks to a calibration pattern and OpenCV (see Fig. 2,right).

We used a straight shaped microcatheter with an outside diameter of 1.7 French (Headway TM 17 regular Microcatheter, Microvention) that was only pushed into the phantom. The microcatheter could automatically be segmented after subtracting a catheter-free image. Refractive stereoscopic reconstruction was thereafter applied followed by an RBF (Radial Basis Functions) smoothing of the curve. The reconstruction error was estimated to be below 0.05 mm on 1000 images of a motionless catheter. A 3D model of the phantom vasculature was segmented from 3DRA data using [START_REF] Yureidini | Robust RANSAC-based blood vessel segmentation[END_REF][START_REF] Yureidini | Local implicit modeling of blood vessels for interactive simulation[END_REF] and manually registered with the images of both cameras. In order to estimate the imprecision of the reconstruction, we measured the average displacement of the catheter tip at the start of the sequence, while no motion was applied (1000 images): a standard deviation of [0.0075, 0.01750.0457] mm was observed along the X-Y-Z axes (Y and Z axes are respectively along the optical axes of the side and top cameras).

Results and Conclusion

We focused our analysis on abrupt motions. Fig. 3 gives an example of a stick and slip transitions: starting from a static position of the catheter tip, a motion amplitude of 10.1 mm was observed on 8 consecutive images, which gives an average speed of 286.86 mm/s. Over this sequence, the instantaneous speeds, as measured in-between two consecutives images, were : [7.06, 123.74, 485.13, 400.04, 388.54, 388.98, 326.56, 27.75] mm/s. A peek of almost 500 mm/s was observed at the start of the sequence.

These first results demonstrate that we have developed a testbed able to validate numerical models and simulations of the mechanical interaction between catheter and blood vessels. In future work, we will use this experimental ground truth to validate and compare existing simulation algorithms. The testbed has been built using a rigid silicon phantom as it allows to simplify the actual environment and puts the emphasis on the tool model and the interaction model between the tool and the blood vessel. Using a soft silicon phantom will be considered in a future work to validate vessel-deformation models but will require a more complex machinery (regular segmentations of the vasculature will be mandatory to capture the deformations of the soft silicon phantom).

Fig. 2 :Fig. 3 :

 23 Fig. 2: (left) Experimental setup, showing the top and side camera, observing the silicon vascular phantom; (middle) image of the phantom by the top camera; (right) calibration pattern used to compute the equation of the refractive plane (inter-bullet distance=2 mm)