Thanh-Tin Nguyen

Maxime Lhuillier

Self-calibration of omnidirectional multi-cameras including synchronization and rolling shutter

Keywords: Bundle adjustment, self-calibration, synchronization, rolling shutter, multi-camera, structure-from-motion

Highlights

• Deal with consumer 360 cameras and spherical cameras without a privileged direction.

• Initialize the time offsets and intrinsic parameters using monocular structure-from-motion.

• Start multi-camera structure-from-motion with central and global shutter assumptions.

• Refine all parameters including time offsets and line delay thanks to a bundle adjustment.

• Experiment on long video sequences using helmet-held multi-cameras.

Introduction

Multi-cameras built by fixing together several consumer cameras become popular thanks to their prices, high resolutions, growing applications including 360 videos (e.g. in YouTube), generation of virtual reality content [1,2], 3D scene modeling [START_REF] Lhuillier | Synchronization and self-calibration for helmetheld consumer cameras, applications to immersive 3d modeling and 360 videos[END_REF]. However such a multi-camera has drawbacks.

First the synchronization of the videos can be a problem. In many cases like GoPro cameras [4], the manufacturer provides a wifi-based synchronization (the user starts all videos at once by a single click). However the resulting time offsets between videos are too inaccurate for applications: about 0.04s and sometimes above 0.1s in our experiments. Assume that a central multi-camera moves at 20km/h (e.g. biking in a city) and two cameras have a time offset equal to only 0.02s, then explain consequences on a 360 video obtained by video stitching. If we neglect this offset, the two videos are stitched as if they have same camera centers at same frame number, although the distance between these centers is 0.11m (20/3.6*0.02). This generates artifacts in the 360 video due to foreground objects that are in the field-of-view (FoV) shared by the two cameras.

Secondly, the low price of a consumer camera implies that the camera is rolling shutter (RS). This means that two different lines of pixels of a frame are acquired at different instants. In a global shutter (GS) camera, all pixels have the same time. If we do a GS approximation of a RS camera, we assume that the camera poses are the same for all lines of a frame although they are not. This can degrade the quality of results in applications such as 3D reconstruction [START_REF] Geyer | Geometric models of rolling-shutter cameras[END_REF], similarly as an inaccurate synchronization degrades the quality of 360 videos.

Last the multi-camera is non-central, i.e. the baseline defined by the distance between the centers of two cameras is not zero. This is inadequate for applications that needs a central multicamera such as 360 video (the smaller the baseline, the better the stitching quality). The user/manufacturer can reduce the baseline (and the multi-camera price) thanks to a small number of cameras. Here we use a DIY multi-camera composed of four Gopro Hero 3 enclosed in a cardboard such that the baseline is as small as possible. Since a small number of cameras also reduces the FoV shared by adjacent cameras, we avoid methods that rely on this shared FoV such as image matching between different monocular videos. A greater number of cameras can be used [6] to increase the shared FoV, but both price and baseline increase. We also experiment using a spherical camera [START_REF]Theta s[END_REF] having only two large FoV (fisheye-like) images.

Our self-calibration takes into account these drawbacks (lack of synchronization, rolling shutter, almost central multicamera) and does not require a calibration pattern. First the multi-camera model is initialized thanks to assumptions that are suitable to an omnidirectional camera without a privileged direction: the cameras have the same setting (frequency, image resolution, FoV) and are roughly equiangular. Second a frame-accurate synchronization is estimated from the instantaneous angular velocities of each camera provided by monocular structure-from-motion. Third both inter-camera poses and intrinsic parameters are refined using multi-camera structurefrom-motion and bundle adjustment. Last we introduce a bundle adjustment that estimates not only the usual parameters but also the sub-frame-accurate synchronization and the rolling shutter. We experiment using videos taken by multi-cameras mounted on a helmet and moving along trajectories of several hundreds of meters or kilometer, then compare our selfcalibration results with ground truth.

Sec. 2 briefly overviews previous work for each step of our method and presents our contributions. Several abbreviations are used in the paper: BA (bundle adjustment), SfM (structure-from-motion), GS (global shutter), RS (rolling shutter), FA (frame-accurate), SFA (sub-frame-accurate), IAV (instantaneous angular velocity), FoV (field-of-view), FpS (frameper-second).

Previous work

Initializing the intrinsic parameters

The intrinsic parameters of a monocular perspective camera can be estimated without a calibration pattern using three steps [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]: projective reconstruction from the given images, selfcalibration assuming that pixels are squares, and refinement using BA [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. If the camera is an axially symmetric fisheye with an approximately known FoV angle, two radial distortion parameters can also be estimated [START_REF] Micusik | Structure from motion with wide circular field of view cameras[END_REF] (this extends the one radial parameter case [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF]). The initialization of these intrinsic parameters is not the paper topic. Here we initialize them assuming that the monocular cameras are roughly equiangular with an approximately known FoV. This is sufficient to experiment our contribution (synchronization and bundle adjustment) and we expect that a method like [START_REF] Micusik | Structure from motion with wide circular field of view cameras[END_REF] improves the results.

Initializing the time offsets

Audio-based synchronization is possible if a distinct sound is available (e.g. a clap) and if the cameras do not have audio/video synchronization issues [1]. A survey of methods for video-based synchronization can be found in [START_REF] Gaspar | Synchronization of two independently moving cameras without feature correspondences[END_REF], but these require inter-camera matching or shared FoV or are designed for non-jointly moving cameras. In our case, we benefit by the assumption of jointly moving cameras but the shared FoV can be too small to automatically obtain a decent matching between two cameras. In [START_REF] Spencer | Temporal synchronization from camera motion[END_REF], transformations are estimated between consecutive frames of every video instead of trying to match different videos. The estimated offset is the one that best "compares" the transformations between two videos. One intuitive example is the translation magnitude that is estimated from tracked features: the larger the translation in one video, the larger the translation in the other. However the transformations in [START_REF] Spencer | Temporal synchronization from camera motion[END_REF] are heuristic (translation) or uncalibrated (homography/fundamental matrix) without radial distortion. Here we propose to compare the IAV estimated by a monocular SfM (Sec. 2.1), which does not have the above inconveniences.

Initializing the inter-camera poses

Once a 3D reconstruction is obtained for every camera (Sec. 2.1) and the time offsets are known (Sec. 2.2), the reconstructions are registered in the same coordinate system. In [START_REF] Carrera | SLAM-based extrinsic calibration of a multi-camera rig[END_REF], a similarity transformation is robustly estimated between two reconstructions using a 3-point RANSAC algorithm and image matching for 3D points in different reconstructions. In [START_REF] Esquivel | Calibration of a multi-camera rig from non-overlapping views[END_REF], the relative pose between two cameras is directly estimated from the pose sequences of their two reconstructions (if the camera motion is not a pure translation). Averaging rotation (e.g. [START_REF] Dai | Rotation averaging with application to camera-rig calibration[END_REF]) can also be used if there is a non constant relative pose between two reconstructions due to the drift of reconstruction(s). The initialization of the inter-camera poses is not the paper topic. Here we initialize them assuming that the multi-camera is roughly central with approximately known inter-camera poses: n cameras that are symmetrically mounted around a symmetry axis. This is enough to feed the bundle adjustment in our cases and [START_REF] Esquivel | Calibration of a multi-camera rig from non-overlapping views[END_REF][START_REF] Carrera | SLAM-based extrinsic calibration of a multi-camera rig[END_REF] can solve this step in all cases.

Global shutter multi-camera bundle adjustments

The BA in [START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF] refines the relative poses between the cameras in addition to the usual parameters (poses of the multi-camera and 3D points) by minimizing a reprojection error. However the reprojection error is in the undistorted space of the classical polynomial distortion model [START_REF] Sturm | Camera models and fundamental concepts used in geometric computer vision[END_REF]. This is due to the fact that the forward-projection of this camera model does not have a closed-form. The BA in [START_REF] Schneider | Bundle adjustment and system calibration with points at infinity for omnidirectional cameras[END_REF] deals with points at infinity, uses ray directions as observations, and transfers the uncertainty from the measure image space to the ray space. The refinement of intrinsic parameters is left as future work in [START_REF] Carrera | SLAM-based extrinsic calibration of a multi-camera rig[END_REF][START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF][START_REF] Schneider | Bundle adjustment and system calibration with points at infinity for omnidirectional cameras[END_REF].

Our multi-camera BA also refines intrinsic parameters (not only inter-camera poses and the other 3D parameters) and minimizes the reprojection error in the right space: the distorted space where the image points are detected. Under the standard assumption that the image noise due to point detection follows zero-mean normalized identical and independent Gaussian vectors, our BA is the Maximum Likelihood Estimator (this assumption is not true in the undistorted space, especially in case of large distortions between undistorted and distorted spaces).

Rolling shutter bundle adjustments

Previous monocular BAs estimate the RS assuming that the 3D points are known in a calibration pattern [START_REF] Oth | Rolling shutter camera calibration[END_REF] or enforce a known RS coefficient [START_REF] Hedborg | Rolling shutter bundle adjustment[END_REF][START_REF] Duchamp | Multiple view 3D reconstruction with rolling shutter cameras[END_REF]. In the context of visual SLAM [START_REF] Klein | Parallel tracking and mapping on a camera phone[END_REF], GS BA is applied to RS (monocular) camera thanks to RS compensation: this method corrects beforehand the RS effects on the feature tracks by estimating instantaneous velocities of the camera. The previous multi-camera BAs estimate neither synchronization nor RS; only [START_REF] Klingner | Street view motion-fromstructure-from-motion[END_REF][START_REF] Saurer | Sparse to dense 3d reconstruction from rolling shutter images[END_REF] deal with known RS but need other sensors.

Every RS BA has a model of the camera trajectory, which provides the camera pose at each instant corresponding to each line of a frame, and which should have a moderated number of parameters to be estimated. In [START_REF] Hedborg | Rolling shutter bundle adjustment[END_REF], one pose is estimated at each frame by BA and the poses between two consecutive frames are interpolated from the poses of these two frames. The BA in [START_REF] Duchamp | Multiple view 3D reconstruction with rolling shutter cameras[END_REF] adds extra parameters to avoid this linear interpolation assumption: it not only optimizes a pose but also rotational and translation speeds at every keyframe. In [START_REF] Oth | Rolling shutter camera calibration[END_REF], a continuous-time trajectory model is used using B-splines and the BA optimizes the knots of the splines. The method chooses the number of knots and initializes their distribution along the trajectory sequence. In [START_REF] Klingner | Street view motion-fromstructure-from-motion[END_REF], the relative pose between an interframe pose and an optimized frame pose is provided by IMU at high frequency. In [START_REF] Saurer | Sparse to dense 3d reconstruction from rolling shutter images[END_REF], rotational and translation speeds are also estimated at every frame (the FpS is only 4Hz) and the BA enforces a relative pose constraint using GPS/INS data. The visual-only RS approaches [START_REF] Klein | Parallel tracking and mapping on a camera phone[END_REF][START_REF] Hedborg | Rolling shutter bundle adjustment[END_REF][START_REF] Duchamp | Multiple view 3D reconstruction with rolling shutter cameras[END_REF][START_REF] Oth | Rolling shutter camera calibration[END_REF] are experimented on few meters long camera trajectories. Our approach is also visual-only and deals with quite longer trajectories (hundreds of meters, kilometers) since it only estimates poses at keyframes.

Self-calibration and synchronization of sensors

In the context of a general multi-sensor, [START_REF] Furgale | Unified temporal and spatial calibration for multi-sensor systems[END_REF] simultaneously estimates the temporal and spatial registrations between sensors. In the experiments, the multi-sensor is composed of a camera and IMU. The best accuracy is obtained thanks to the use of all measurements at once, a continuous-time representation (a B-spline for IMU poses) and maximum likelihood estimation of the parameters (time offset, transformation between IMU and camera, IMU poses, and others). In [START_REF] Lovegrove | Spline fusion: a continuoustime representation for visual-intertial fusion with application to rolling shutter cameras[END_REF], a camera-inertial multi-sensor is self-calibrated (synchronization, spatial registration, intrinsic parameters) by a sliding window visual odometry. Thanks to an adequate continuous-time motion parametrization, it also deals with RS cameras and has a better parametrization of the rotations. Indeed, it avoids the singularities of the global and minimal parametrization of rotations (e.g. in [START_REF] Furgale | Unified temporal and spatial calibration for multi-sensor systems[END_REF]), but assumes that the time between consecutive keyframes is uniform. Our work introduces a global minimal rotation parametrization and deals with non-uniform distribution of keyframes provided by standard SfM [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF].

Recently, [START_REF] Vo | Spatiotemporal bundle adjustment for dynamic 3d reconstruction[END_REF] synchronizes and self-calibrates consumer cameras using BA in a different context: assumptions are removed (rigidity on both multi-camera and scene), others are added (FoV shared by cameras, physics-based motion priors for moving objects), and the rolling shutter is not estimated.

Our contributions

Our multi-camera BA estimates the SFA synchronization and the line delay coefficient of the RS. Furthermore, this is done over long video datasets without additional sensors (hundreds of meters or kilometers). The previous work do not do this. In contrast to [START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF], our BA also estimates the intrinsic parameters and minimizes the reprojection errors in the original image space (not the rectified one) with the same polynomial distortion model [START_REF] Sturm | Camera models and fundamental concepts used in geometric computer vision[END_REF][START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration ?[END_REF]. Another contribution is the FA synchronization that deals with cameras with small/empty shared FoV. As mentioned in Secs. 2.1 and Sec. 2.3, our initialization does not intend to compete with the accuracy and generality of previous initializations of intrinsic parameters and relative poses.

Contributions over our previous conference work [START_REF] Lhuillier | Synchronization and self-calibration for helmetheld consumer cameras, applications to immersive 3d modeling and 360 videos[END_REF][START_REF] Nguyen | Adding synchronization and rolling shutter in multi-camera bundle adjustment[END_REF] are the following: check approximations [START_REF] Klingner | Street view motion-fromstructure-from-motion[END_REF] in the computation of the image projection that takes into account RS and synchronization, refine simultaneously calibration and synchronization/RS, deal with spherical cameras, more details on synchronization and BA (rotation parametrization, sparsity of solved system, derivatives of implicit reprojections errors). There are also new experiments on synchronization (comparison with ground truth, robustness of SFA refinement with respect to bad FA initialization), stability of both SFA synchronization and RS over time in long videos and with respect to keyframe sampling.

Overview of our algorithm

First the monocular camera model (we experiment the classical polynomial distortion model [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration ?[END_REF][START_REF] Sturm | Camera models and fundamental concepts used in geometric computer vision[END_REF][START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF] and the unified camera model [START_REF] Geyer | A unifying theory for central panoramic systems and practical implications[END_REF]) is initialized in Sec. 4 assuming that the fisheyes are roughly equiangular and using an approximate knowledge of their FoV angle.

Second we apply monocular SfM [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] and calibration refinement by BA for every camera. However SfM can fail for a video due to the combination of two difficulties: lack of texture and approximate calibration. We assume that there is at least one textured enough video such that this is successful. Since the cameras have the same setting, we benefit by the refined intrinsic parameters by BA to redo the monocular SfM of the other videos. Thus a difficulty (approximate calibration) is reduced for the less textured videos and the risk of failure decreases.

Third the FA synchronization between all videos is obtained by using the method in Sec. 5. We skip few frames in each video such that the sequels of the videos are FA synchronized: from now frames with the same index are taken at the same time up to the inverse of the FpS.

Fourth a central multi-camera calibration is initialized from the estimated intrinsic monocular parameters and approximate inter-camera rotations (Sec. 2.3).

Fifth we apply multi-camera SfM [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] followed by multicamera BA [START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF] by adding the intrinsic parameters as new estimated parameters. Up to now, we did three approximations: global shutter, central multi-camera, and zero sub-frame residual time offsets. Furthermore we only applied the SfMs (both monocular and multi-cameras) on the beginning of the videos to obtain initial synchronization and calibration (the 2k first frames in our experiments). Then the multi-camera SfM is applied a second time on the whole videos.

Last we apply the multi-camera BA in Sec. 6 for estimating the SFA synchronization and the line delay with usual parameters.

Sec. 7 explains how to efficiently compute non-closed form image projections and their derivatives involved in BA. The experiments and conclusion are in Secs. 8 and 9, respectively.

Equiangular initializations

Secs. 4.1 and 4.2 describe two monocular camera models and their initializations (before all SfM and BA computations). Both models involve the intrinsic parameter matrix K of a perspective camera: K has focal parameters f x and f y , principal point z 0 and zero skew. The classical polynomial distortion model [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration ?[END_REF][START_REF] Sturm | Camera models and fundamental concepts used in geometric computer vision[END_REF][START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF] is often used since its closed-form backprojection is useful for SfM tasks and epipolar geometry. It has several radial distortion parameters and can be applied to consumer cameras like Gopro [4]. The unified camera model [START_REF] Geyer | A unifying theory for central panoramic systems and practical implications[END_REF] is also interesting since it deals with fisheyes having FoV larger than 180 • (like those of spherical camera [START_REF]Theta s[END_REF]) although it only has a single radial distortion parameter. The equiangular initialization can be adapted to other camera models.

Polynomial distortion model

Back-projection

The function from the distorted (i.e. original) image to the undistorted (i.e. rectified) image depends on radial distortion parameters k i (tangential distortions are neglected). Let z d and z u be the distorted and undistorted coordinates of a pixel. Their normalized coordinates zd and zu meet

K zd 1 = z d 1 and K zu 1 = z u 1 . (1)
Let rd = ||z d || be the normalized radial distance in the distorted image. The relation between distorted and undistorted coordinates is

zu = (1 + n i=1 k i r2i d )z d . (2)
Lastly, the back-projected ray of pixel z d has direction z u 1 in the camera coordinate system.

Initialization

Here we initialize k i , z 0 , f x and f y for an equiangular camera. The camera is equiangular if the angle µ between the principal direction 0 0 1 and the back-projected ray is proportional to the (non-normalized) radial distance r d in the distorted image. We have

r d = ||z d -z 0 || and tan µ = ||z u ||. If the camera is equiangular, f x = f y = f and there is a constant c such that µ = cr d . Thus µ = c f rd . Since ||z u || = rd (1 + n i=1 k i r2i d ), tan(c f rd ) = tan µ = rd + n i=1 k i r2i+1 d . ( 3
)
Since tan is not a polynomial, Eq. 3 can not be exact. We use a Taylor's approximation

tan µ ≈ n i=0 t i µ 2i+1 = µ + µ 3 3 + 2µ 5 5 + 17µ 7 315 + • • • (4)
and identify coefficients between Eqs 3 and 4. We obtain c f = 1 using t 0 and k i = t i if i ≥ 0. In practice, we initialize z 0 at the image center and compute f = r d /µ for a pixel z d at the center of an image border where the half-FoV µ is approximately known.

Unified camera model

Forward projection

Let x = x y z ∈ R 3 \ {0} be a 3D point in the camera coordinate system. Let S be the unit sphere in R 3 centered at 0 and let ξ ∈ R + . The projection p(x) of x by this model is obtained as follows: first x is projected onto S, then x/||x|| is projected onto the image plane by a perspective camera with the center 0 0 -ξ and the intrinsic parameter matrix K. Formerly,

p(x) = π(K( x ||x|| +           0 0 ξ          
)) where π(

          u v w           ) = u/w v/w . (5)

Initialization

Here we initialize ξ, z 0 , f x and f y for an equiangular camera. Let µ be the angle between the principal direction 0 0 1 and the back-projected ray, which is a half-line started at 0 with the direction x/||x||. Appendix A shows that

f x = f y = f ⇒ ||p(x) -z 0 ||/ f = sin µ ξ + cos µ . ( 6
)
If the camera is equiangular, f x = f y = f and there is a constant c such that µ = c||p(x)z 0 ||. Since sin µ ξ+cos µ is not linear in µ, we approximate it thanks to Taylor's expansions:

sin µ ξ + cos µ ≈ µ -µ 3 /6 ξ + 1 -µ 2 /2 = µ(1 -µ 2 /6) (1 + ξ)(1 -µ 2 /(2ξ + 2)) = µ 1 + ξ (1 + µ 2 ( 1 2ξ + 2 - 1 6 ) + O(µ 4 )). ( 7
)
We initialize ξ = 2 such that this approximation is linear in µ. Now we distinguish two cases for the initialization of z 0 and f . If every pixel of the (rectangular) image has a back-projected ray, we initialize z 0 at the image center and take a point z 1 at the center of an image border where the half-FoV µ is approximately known. Otherwise, we assume that the pixels that have back-projected rays form a disk whose radius and center can be estimated. Then we initialize z 0 by this center and take a point z 1 at the disk boundary where the half-FoV µ is approximately known. In both cases, f is initialized by Eq. 6 using p(x) = z 1 .

Synchronization initialization

The synchronization initialization is required by the multicamera SfM-BA and has two steps. First Sec. 5.1 estimates instantaneous angular velocities thanks to monocular SfM-BA and global shutter approximation. Then time offsets are computed by correlation of IAVs of different cameras; Secs. 5.2 and 5.3 describe the two-and multi-camera cases respectively.

Instantaneous angular velocity (IAV)

Every monocular video is reconstructed such that every frame has a computed pose (both keyframes and nonkeyframes, not only keyframes as in the paper remainder). Thus the keyframe-based SfM [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] is followed by pose calculations for the non-keyframes and by BA. In practice, it is sufficient to reconstruct few thousands of frames at the video beginning for the synchronization initialization. Let R t i be the rotation of the pose of the t-th frame in the reconstruction of the i-th video. The IAV θ t i at the t-th frame (of the i-th video) is approximated by the angle of rotation R t+1 i (R t i ) , i.e.

θ t i = arccos((trace(R t+1 i (R t i ) ) -1)/2). ( 8
)
We omit the FpS coefficient since all cameras have the same.

Intuitively, two frames of different but jointly moving cameras have same IAV if they are taken at the same time. This is shown in Appendix B by taking account the fact that the R t i are expressed in arbitrary coordinate systems due to the monocular SfM.

Synchronize two cameras

We compute an IAV table for every camera and find the time offset that maximizes the correlation (ZNCC) between two such tables (match two sub-tables with the same length in different tables). The time offset o i, j between the i-th and j-th cameras maximizes correlation ZNCC i, j between vectors θ t i and θ t+o i, j j . We also introduce a simple SFA refinement method. The sub-frame offsets are estimated like sub-pixelic disparity using a quadratic fit [START_REF] Szeliski | Symmetric sub-pixelic stereo matching[END_REF]: first approximate the function from o i, j to ZNCC i, j using a quadratic polynomial defined by its 3 values at o i, j + {-1, 0, +1}; then estimate i, j such that o i, j + i, j maximizes this polynomial. In contrast to the FA offsets o i, j ∈ Z, the SFA offsets o i, j + i, j ∈ R are not used for the input of our BA.

Consistently synchronize more than two cameras

We remind that the goal of the FA synchronization is to skip s i frames at the beginning of the i-th video such that the sequels of the videos are FA synchronized (this is required for multicamera SfM). Thus o i, j = s js i for all i j, which in turn imply that the sum of the offsets along every loop in the camera graph should be zero (e.g. we should have o 0,1 + o 1,2 + o 2,0 = 0 for loop 0 → 1 → 2 → 0). However such a sum can be nonzero since the offsets are estimated independently.

There are several ways to deal with this loop constraint. First only compute offsets o 0,i . But this solution privileges a camera. Second compute all offsets o i, j , generate candidate offsets around o i, j for every pair (i, j), and select the candidate offsets that maximizes i j ZNCC i, j such that the sum of candidate offsets along every loop is zero. We implement an intermediate and simple solution where every camera has the same importance assuming that the cameras are symmetrically mounted around a symmetry axis: we only consider the spatial adjacency of the n cameras, i.e. we only compute offsets

o 0,1 , o 1,2 , • • • o n-2,n-1 , o n-1,0 (instead of all o i, j ) and only use the loop 0 → 1 → • • • n -1 → 0 (instead
of all loops) in the scheme above. In practice, we found that it is sufficient to generate candidate offsets that differ from the initial ones by +1 or 0 or -1.

In the remainder of the paper, we use the notation o i, j for offsets that meet the loop constraint.

Bundle adjustment for RS and synchronization

This is the last step of our method and it requires the multicamera initialization described in Sec. 6.1. Sec. 6.2 presents our continuous-time parametrization of the multi-camera motion: it is defined by the composition of a function M from a time interval to R 3 × R k and a function R from R k to the set of rotations in R 3 . Sec. 6.3 describes a keyframe of the multi-camera, where every line has a time that depends on the camera that captures the line, its y-coordinate and the line delay. Sec. 6.4 approximates M(t) at a time t from the few M(t i ) corresponding to the beginnings of the keyframes; this is useful to moderate the number of parameters estimated by BA. Sec. 6.5 provides a simple method to compute the reprojection error minimized by BA. Last Secs. 6.6 and 6.7 are more technical: we choose R in the former and detail the sparse structure of the linear system solved by BA in the latter.

Initialization

First we assume that the monocular videos are FA synchronized by removing few frames at their beginning (Sec. 5). Then we define the i-th frame of the multi-camera by a concatenation of sub-images, every of them is the i-th frame of a monocular camera. From now on, we use word frame for "frame of the multi-camera" and the video is the sequence defined by all these frames. Last we use a standard SfM based on keyframe sub-sampling of the video (Appendix G) and local BA [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] assuming GS. We remind that the keyframes are the only frames whose poses are refined by the BAs (this is useful for both time computation and accuracy).

Parametrization of the multi-camera trajectory

Let R be a C 1 continuous and surjective function that maps R k to the set of the 3D rotations (typical values are k ∈ {3, 4}). We assume that there is a C 3 continuous function M : R → R 3 ×R k that parametrizes the motion of the multi-camera. More precisely, M(t) T = T M (t) T E M (t) T where t ∈ R is the time, T M (t) ∈ R 3 is the translation and R(E M (t)) is the rotation of the multi-camera pose. The columns of R(E M (t)) and T M (t) are the vectors of the multi-camera coordinate system expressed in world coordinates. The choice of R (including E M and k) is detailed in Sec. 6.6 for the paper clarity.

Thanks to these notations and assumptions, we will approximate M(t) by using values of M taken at few times t 1 , • • • , t m . Then our model of the camera trajectory not only provides the multi-camera pose at each instant corresponding to each line of a frame, but it also has a moderated number of parameters to be estimated by BA: the vector concatenating all M(t i ), which has dimension m(3 + k). Sec. 6.3 defines t i and Sec. 6.4 describes our approximations of M(t) by using the M(t i ).

Time, RS and synchronization parameters

The i-th keyframe is an image composed of sub-images taken by the monocular cameras. Every line of every sub-image is taken at its own time, which is described now. The 0-th line of the 0-th sub-image in the i-th keyframe is taken at time t i , assuming that the time exposure of a line is instantaneous [START_REF] Geyer | Geometric models of rolling-shutter cameras[END_REF]. Thus t i+1t i is a multiple of the inverse of the FpS. Since the cameras are RS, the line delay τ is such that the y-th line of the 0-th sub-image in the i-th keyframe is taken at time t i + yτ. Let ∆ j ∈ R be the sub-frame residual time offset between the j-th video and the 0-th video. Then the 0-th line of the j-th sub-image in the i-th keyframe is taken at time t i + ∆ j . Since we assume that all cameras have the same FpS and same (and constant) τ, the y-th line of the j-th sub-image in the i-th keyframe is taken at time t i + ∆ j + yτ. Fig. 1 illustrates the trajectory M(t) of a multi-camera defined by four monocular rolling shutter cameras having non-zero time offsets ∆ j .

Approximations for the multi-camera trajectory

Let ∆ = max i (t i+1t i ) and shortened notation m i = M(t i ). Thanks to the C 3 continuity of M and Taylor's expansions of M at t i , we explicit two approximations M 1 (t) and M 2 (t) of M(t) in the neighborhood of t i as functions of m i-1 , m i and m i+1 . These approximations have remainders expressed in terms of ∆ and |tt i |. By neglecting these remainders, we compute M(t) for the y-th line of the j-th camera/sub-image in the i-th keyframe using t = t i + ∆ j + yτ (Sec. 6.3) during our BA.

Linear approximation M 1 of M

We have Taylor's linear expansion

M(t) = m i + (t -t i )M (t i ) + O(|t -t i | 2 ) (9)
and express the derivative M (t i ) as a function of m i-1 , m i and m i+1 . Let reals a > 0 and b > 0, vectors x, y, z in R k+3 , function

D 1 such that D 1 (x, y, z, a, b) = bz a(a + b) - ax b(a + b) + (a -b)y ab , (10)
and shortened notation

D i 1 = D 1 (m i-1 , m i , m i+1 , t i+1 -t i , t i -t i-1 ). (11)
Appendix C shows that

M (t i ) = D i 1 + O(∆ 2 )
. We obtain

M 1 (t) = m i + (t -t i )D i 1 if t ≈ t i . ( 12
) If i = 0 (similarly if i = m), we use D 0 1 = m 1 -m 0 t 1 -t 0 .
6.4.2. Quadratic approximation M 2 of M Similarly, we have Taylor's quadratic expansion of M at t i and express the derivative M (t i ) as a function of m i-1 , m i and m i+1 . Let a, b, x, y, z as in Sec. 6.4.1, function D 2 such that

D 2 (x, y, z, a, b) = 2z a(a + b) + 2x b(a + b) - 2y ab , (13)
and shortened notation

D i 2 = D 2 (m i-1 , m i , m i+1 , t i+1 -t i , t i -t i-1 ). ( 14
) Appendix C shows that M (t i ) = D i 2 + O(∆). We obtain M 2 (t) = m i + (t -t i )D i 1 + (t -t i ) 2 2 D i 2 if t ≈ t i . ( 15
) If i = 0 (similarly if i = m), we use D 0 1 = m 1 -m 0 t 1 -t 0 and D 0 2 = 0.

Reprojection error of the multi-camera

Since our BA minimizes the sum of squared modulus of reprojection error for every inlier, this section describes the computation of a reprojection error for a 3D point x ∈ R 3 (in world coordinates) and its inlier observation p ∈ R 2 in the j-th subimage of the i-th keyframe.

First we introduce notations. Let p ∈ R 2 be the projection of x in the j-th sub-image of the i-th keyframe. The reprojection error is p -p. Let (R j , t j ) be the pose of the j-th camera in the multi-camera frame. Let p j : R 3 \ {0} → R 2 be the projection function of the j-th camera. We assume that p j , R j , t j are constant. The acquisition times of p = (x, y) and p = ( x, ỹ) are

t p = t i + ∆ j + yτ and t p = t i + ∆ j + ỹτ. ( 16
)
Second we detail the relation between p and x. Both E M (t p ) and T M (t p ) (i.e. M(t p ) in Sec. 6.2) are defined by one equation chosen among Eq. 12 and Eq. 15 using the index i of the keyframe and t = t p . The coordinates of x in the multi-camera coordinate system is

x M = R(E M (t p )) (x -T M (t p )). ( 17
)
The coordinates of x in the j-th camera coordinate system and the projection of x are

x j = R j (x M -t j ) and p = p j (x j ). ( 18
)
We see that p needs the computation of x M , which in turn needs the computation of (the y coordinate of) p. Such a problem is solved thanks to an approximation in [START_REF] Klingner | Street view motion-fromstructure-from-motion[END_REF]: t p is replaced by t p in Eq. 17, i.e. we assume that the multicamera pose is the same at times t p and t p . We think that this is acceptable since |t pt p | ≤ τ||p -p|| and the magnitude order of τ is 10 -5 s/pixel and p is an inlier (i.e. ||p -p|| ≤ 4 pixels). Sec. 7.3 presents another solution without this approximation.

Parametrization of rotations

Sec. 6.6.1 lists useful properties of R (reminder: R is a function introduced in Sec. 6.2 which maps R k to the set of the 3D rotations). Then Sec. 6.6.2 explains our choice of R to meet the properties in Sec. 6.6.1.

Details on R properties

First we note that R is a global parametrization used for the whole camera trajectory (we do not use local parametrizations, i.e. different parametrizations for different keyframes). Second the C 1 continuity of R is needed by BA for the derivative computations of the reprojection errors. Third we follow [START_REF] Oth | Rolling shutter camera calibration[END_REF] by using a minimal (non-redundant) parametrization R of the rotations to limit the number of estimated parameters. Thus k = 3.

Fourth BA needs another property. According to Secs. 6.2 and 6.3,

m i = M(t i ) = T M (t i ) E M (t i ) ∈ R 6 is one of the parameter vectors estimated by BA such that (T M (t i ), R(E M (t i )))
is the pose (of the first line) of the i-th keyframe. Since the set of all rotations in a neighborhood of a current estimate of rotation R(E M (t i )) should be reachable by the parametrization R during every BA iteration [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF], the jacobian ∂R of R should be rank 3 at E M (t i ). In other words, E M (t i ) should not be a singularity of R.

Unfortunately, every 3D parametrization R of the rotation set has singularities [START_REF] Singla | How to avoid singularity when using Euler angles ?[END_REF]. Thus we choose R in Sec. 6.6.2 such that all its singularities are far from the multi-camera motion that we want to refine using BA.

Choice of a minimal parametrization R

First we consider R candidates and describe constraints that they induce on a class of multi-camera motions: all yaw motions are possible but pitch and roll are small. Such motions are very common for a helmet-held multi-camera and an user exploring the environment without special objective like grasping at object on the ground (and also for a car-fixed multi-camera).

Even the popular exponential map ω → exp([ω] × ) has singularities: they form concentric spheres with center 0 and radii that are multiples of 2π [START_REF] Grassia | Practical parametrization of rotations using the exponential map[END_REF]. This can be seen thanks to the equivalent angle-axis (θ, n) representation where ||n|| = 1 and ω = θn. Thus the range of the angle θ is equal to 4π for every axis n. If we choose R(ω) = exp([ω] × ) as in [START_REF] Oth | Rolling shutter camera calibration[END_REF] and would like to avoid the singularities, the multi-camera should avoid multiple turns on the left (or right) around buildings and avoid straight trajectory segments where ||ω|| ≈ 2πk and k ∈ Z * .

We also detail the case of Euler's parametrization

E(α, β, γ) = R z (γ)R y (β)R x (α) (19)
where R x (α), R y (β) and R z (γ) are the rotations about respective axes 1 0 0 , 0 1 0 , 0 0 1 and with respective angles α, β, γ. The singularities (α, β, γ) of E form parallel and equidistant planes of equations β = π/2 + pπ such that p ∈ Z [START_REF] Singla | How to avoid singularity when using Euler angles ?[END_REF]. If we choose R = E and the coordinate systems (both world and multi-camera) are such that ∀i, R(E M (t i )) ≈ R x (α i ), we are far from the singularities. If the coordinate systems are such that ∀i, R(E M (t i )) ≈ R y (π/2)R x (α i ), we are close to the singularities.

Last we choose R inspired by the Euler's case above. Let

R(α, β, γ) = AR z (γ)R y (β)R x (α)B ( 20
)
where rotations A and B do not depend on (α, β, γ). We estimate A and B such that β is close to 0 for all keyframe rotations of the multi-camera trajectory before the BA in Sec. 6 (technical details in Appendix D). Now the camera motion in our class is far from all singularities. Note that the local Euler parametrization, that is used in BA [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF], is a special case of this parametrization.

Sparsity of the reduced camera system (RCS)

In this section, we explicit the sparsity of the RCS that is solved by BA [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. This is important for efficient computations.

Notations and global structure of the RCS

The m vectors m i = M(t i ) are the parameters of the multicamera trajectory (Sec. 6.2), they meet m i ∈ R 6 (Sec. 6.6.1), and we define

M = m 1 • • • m m ∈ R 6m . Let m ∈ R m
be the other optimized camera parameters among intrinsic parameters, camera poses in multi-camera coordinates, line delay and time offsets. Since these other optimized parameters are the sames at all keyframes, m 6m. For example, m = 1000 and m ≤ 4+4 * 15 = 64 if the multi-camera has four Gopro cameras: there are line delay τ, time offsets ∆ 1 , ∆ 2 , ∆ 3 (Sec. 6.3), and every camera has parameters f x , f y , u 0 , v 0 , k 1 • • • k 5 (Sec. 4.1.1) and 6D pose in multi-camera coordinates. Let X be the vector that concatenates points x l ∈ R 3 in world coordinates. The projection function of x l in the i-th keyframe is concisely written ϕ(m i-1 , m i , m i+1 , m , x l ) for both approximations in Sec. 6.4 (omit m i-

1 if i = 1 and omit m i+1 if i = m). Let H = ∂ϕ ∂M ∂ϕ ∂m ∂ϕ ∂X ∂ϕ ∂M ∂ϕ ∂m ∂ϕ ∂X =           U U W (U ) U W W (W ) V           . ( 21
)
Here H is the approximated hessian of the cost function minimized by BA. It is defined as a sum for all 2D inliers (detailed notations are omitted). The RCS is

U U (U ) U - W W V -1 W W = Z Z (Z ) Z . (22)
We have Z = U -WV -1 W ∈ R 6m×6m , Z ∈ R 6m×m and Z ∈ R m ×m . Since m 6m, Z is the preponderant block in the RCS and we only focus on the Z sparsity.

Sparsity of Z

Here we represent Z by a shape included in Z 2 , i.e. a set of pixels in an image such that every pixel corresponds to a nonzero 6 × 6-block of Z. Then we show in Appendix E that the shape of our Z (which involves SFA and RS using projection function ϕ(m i-1 , m i , m i+1 , m , x l )) is included in a dilation of the shape of the standard Z (which involves FA and GS using projection function ϕ(m i , m , x l )) by {-1, 0, +1} 2 . We remind that this dilation is an operation morphology that expands a shape by one pixel in both dimensions and both directions. Thus our RCS is slightly less sparse that the standard RCS (in practice it is very similar according to the example in Fig. 2).

In the case where the loops are not closed in the video and the track length is bounded by l, the standard Z is a 6 × 6-blockwise band matrix with bandwidth l (Sec. A6.7.1 in [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]) and our Z is a 6 × 6-block-wise band matrix with bandwidth l + 1.

Non-closed-form image projections

Sec. 7 computes the projection p and its derivatives with respect to a vector θ of parameters optimized by BA, although p does not have a closed-form expression from θ. This is needed by BA and occurs in two cases in this paper. The general case is described in Sec. 7.1 assuming that p and θ meet an equation that implicitly defines p for a given value of θ. Then Secs. 7.2 and 7.3 apply the general case in two cases: the projection of the polynomial distortion model (reminder: Sec. 4.1.1 only computes the back-projection) and the exact projection using a continuous-time trajectory model (reminder: Sec. 6.5 only computes an approximate projection). This problem and its solution are similar to those of a general non-central catadioptric camera (Appendix D in [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF]).

General case

We know an approximate value p of p ( p is an inlier detected in an image), a C 1 continuous function g(z, θ) from R 2 × R p to R 2 such that p is the solution z of g(z, θ) = 0, and the current value θ 0 of θ (provided by initialization or previous iteration of BA). First p is estimated by non-linear least-squares minimizing z → ||g(z, θ 0 )|| 2 . In practice, we use the iterative Gauss-Newton's method starting from z = p with no more than 5 iterations (Newton's method can also be used). Then the implicit function Theorem implies that we locally have a C 1 continuous function ψ such that p = ψ(θ) if det ∂g ∂z 0. By differentiating g(ψ(θ), θ) = 0 using the Chain rule, we obtain

∂p ∂θ = ∂ψ ∂θ = -( ∂g ∂z ) -1 ∂g ∂θ . (23)

Case 1: polynomial distortion model

Here we focus on the projection p = p j (x j ) in Eq. 18 using the camera model in Sec. 4.1.1 and assuming that the 3D point x j is known (in camera coordinates). First we define θ and g by

θ = ( f x , f y , z 0 , k 1 , k 2 , • • • , k n , x j ), z u = π(Kx j ), (24)
u v = z, u 0 v 0 = z 0 , r2 = (u -u 0 ) 2 f 2 x + (v -v 0 ) 2 f 2 y , (25)
g(z, θ) = (1

+ n i=1 k i r2i )(z -z 0 ) -z u + z 0 . (26)
Then we show that g(p, θ) = 0. Since z u 1 in Sec. 4.1.1 and x j are colinear, z u is the same in Sec. 4.1.1 and Eq. 24.

Furthermore, p = z = z d implies that rd (in Sec. 4.1.1) and r are the same. We obtain g(p, θ) = 0 by multiplying Eq. 2 on the left by f x 0 0 f y . Last we apply Sec. 7.1: find p by minimizing z → ||g(z, θ 0 )|| 2 and use Eq. 23 for p derivatives:

∂p ∂k i = -( ∂g ∂z ) -1 r2i (z -z 0 ), ( 27
)
∂p ∂x j = ( ∂g ∂z ) -1 ∂z u ∂x j . ( 28
)
Using a similar function g, Appendix F shows that

∂p ∂ f x ∂p ∂ f y ∂p ∂u 0 ∂p ∂v 0 =       u-u 0 f x 0 1 0 0 v-v 0 f y 0 1       . (29)
Eqs. 27, 28 and 29 need z = p and θ = θ 0 . We note that the derivative computations in Eq. 28 are easy from those of a standard perspective camera (i.e. ∂z u ∂x j ): multiply on the left by ( ∂g ∂z ) -1 . This also holds for derivatives with respect to the parameters defining x j (Eq. 18) thanks to the Chain rule.

Case 2: exact calculation for RS and synchronization

Here we focus on the projection p by the j-th camera in the i-th keyframe of the 3D point x in world coordinates without the approximation in Sec. 6.5. Let p(θ j , m, x) be the projection of x where m = M(t) is the parameter of the multi-camera pose (reminder: M(t) is introduced in Sec. 6.2) and the vector θ j concatenates the intrinsic/distortion parameters and the pose of the j-th camera in the multi-camera coordinate system. The acquisition time of p is t(∆ j , τ, p) = t i + ∆ j + τ 0 1 p.

(

) 30
and we use notation M(m i-1 , m i , m i+1 , t) for the chosen approximation (Eq. 12 or Eq. 15). Thus we have

p = p(θ j , M(m i-1 , m i , m i+1 , t(∆ j , τ, p)), x). ( 31
)
Now we define θ and g by

θ = (θ j , m i-1 , m i , m i+1 , ∆ j , τ, x), (32)
g(z, θ) = p(θ j , M(m i-1 , m i , m i+1 , t(∆ j , τ, z)), x) -z.( 33
)
We see that g(p, θ) = 0. Then we apply Sec. 7.1: find p by minimizing z → ||g(z, θ 0 )|| 2 and use Eq. 23 for the p derivatives.

If we use the linear trajectory approximation M 1 (Eq. 12), we have a simple expression

∂g ∂z = ∂p ∂m ∂M ∂t ∂t(∆ j , τ, z) ∂z -I 2 (34) = τ ∂p ∂m D i 1 0 1 -I 2 ∈ R 2×2 . ( 35
)
We note that the derivative computations without approximation (Eq. 23) can be deduced from those with approximation (i.e. ∂p ∂θ = ∂g ∂θ ( p, θ 0 )): replace p by p in the derivative by θ and multiply it on the left by -( ∂g ∂z ) -1 .

Experiments

Datasets

Secs. 8.1.1 and 8.1.2 present cameras and video sequences that are used in the experiments, respectively. Tab. 1 summarizes our dataset (both cameras and videos).

Cameras

The consumer multi-cameras are modeled by several rigidly mounted monocular cameras and the user fixes them on a helmet. We assume that all calibrations parameters (time offsets ∆ j , line delay τ, intrinsics, radial distortion, relative poses) are constant during a video acquisition. The camera gain is not fixed and evolves independently for every camera.

First there are 360 cameras composed of four GoPro Hero 3 cameras [4], that are started by a single click on a wifi remote. The user can choose the relative poses of the cameras: they are enclosed in a cardboard for small baseline or are fixed by using the housings provided with the cameras (larger baseline).

Second there is a spherical camera modeled by two opposite fisheyes (no relative pose choice) that are synchronized. The Ricoh Theta S multi-camera [START_REF]Theta s[END_REF] has a very small baseline thanks to the use of a prism mirror in front of a monocular camera (its FoV is split in two equal parts, each of them sees more than a half-sphere as a real fisheye does).

We also experiment on a professional multi-camera (Point-Grey Ladybug 2 [START_REF]Ladybug2[END_REF]) since its ground truth is provided by the manufacturer (as a table of rays) and also for experimenting on an ideal multi-camera with global shutter and perfect synchronization. Except in a synthetic case, the other cameras have incomplete ground truth (a strobe always provides τ).

Videos

There are three real multi-camera videos taken under various conditions using four GoPro cameras (BC1: bike riding in a city, WT: walking in a town, FH: paragliding flying at very low height above a hill). WT is taken in the early morning during summer to avoid moving car and people, but the lighting is low. BC1 has sunny lighting (with contre-jours) and most cars are parked. WT and BC1 have the same calibration setting, FH has larger baseline and lower FpS and better angular resolution. Fig. 3 shows images and cameras of BC1 and WT, Fig. 4 shows images and cameras of FH.

BikeCity2 (BC2) is generated by ray-tracing of a synthetic urban scene having real textures and by moving the camera Name (short name) along a trajectory that mimics that of BC1 (the "pose noise", i.e. the relative pose between consecutive frames, are similar in both videos). We obtain a video for each camera by compressing the output images using ffmpeg and options "-c:v libx264 -preset slow -crf 18". BC2 has ground truth: f ∆ 1 = 0.25, f ∆ 2 = 0.5, f ∆ 3 = 0.75 and similar τ as BC1 (reminder: if f ∆ j = 1 and f is the FpS, ∆ j is the time between two consecutive frames). Fig. 3

Camera f r (mr) b (cm) τ (µs) f ∆ j l (m) f r FoV k f r #Tracks ||β i || ∞ BikeCity1 (BC1) 4*Gopro

also shows images of BC2.

There is also video WU (walking in the campus of the UCA university) using a spherical camera: the Ricoh Theta S (Fig. 5).

Last CarCity is taken by the Ladybug 2 and has a similar trajectory as BC1. However, it is mounted on a car (using a mast) and is about 4 meters above the ground, as shown in Fig. 6. The images are uncompressed (all others are videos compressed using H.264).

Main notations

We use shortened notations:

• #2D= number of 2D inliers

• GT= ground truth,

• f = FpS

• RMS= RMS error in pixels in the original (distorted) image space,

• BA= bundle adjustment.

• method gs.sfa is the SFA refinement in Sec. 5.2

• y max is the number of lines of a monocular image (y max = 768 for CC, y max = 1440 for FH, y max = 960 for BC1,WT,BC2 and WU).

Our BA is named by a combination of several notations that describes the estimated parameters:

• C (central approx.) estimates all rotations R j and fixes all translations t j = 0 (reminder: (R j , t j ) is the pose of the j-th camera in the multi-camera frame)

• NC (non-central) estimates all (R j , t j )

• RS (rolling shutter) estimates the line delay τ

• GS (global shutter) fixes τ = 0

• SFA (sub-frame accurate) estimates all time offsets ∆ j

• FA (frame accurate) fixes all ∆ j = 0

• INT (intrinsics) estimates all intrinsic parameters:

( f x , f y , u 0 , v 0 , ξ) or ( f x , f y , u 0 , v 0 , k 1 • • • k 5 )
depending on the camera model chosen in Sec. 4.2.1 or Sec. 4.1.1 respectively; every camera has its own parameters.

Thus GS.NC.SFA.INT (or gs.nc.sfa.int) is a BA that fixes τ = 0 and estimates simultaneously all ∆ j , R j , t j and intrinsic parameters and keyframe poses m i and 3D points. The threshold for the inlier selection is set to 4 pixels in all videos. Every BA has three inlier updates, each one is followed by a Levenberg-Marquardt minimization for these inliers. A succession of two BAs is possible, e.g. gs.c.fa.int+rs.c.sfa.

The error e(∆) is the sum of the absolute errors of all f ∆ j . The error e(τ) is the relative error of τ. We also define the error of the estimated multi-camera calibration by a single number d, which is the RMS for all multi-camera pixels of the angle between rays of the two calibrations (the estimated one and the GT one) that back-project the same pixel. There are two reasons to do this. First the accuracy is only needed for the ray directions in applications (SfM, video stitching, 3d modeling of a scene) in the central case. Second parameters can compensate themselves if their estimations are biased (e.g. rotation/principal point near-ambiguity for one view [START_REF] Agapito | Self-calibration of rotating and zooming cameras[END_REF]). Now we detail the computation of d. Since the rays of the estimated calibration and the rays of the GT calibration can be expressed in different coordinate systems, we estimate a registration between both coordinate systems before computing angles between rays. The registration is defined by a rotation R, that maps one ray set to the other ray set (ignoring translation of ray origins). More precisely, R is the minimizer of e(R) = N i=1 ||r gt i -Rr est i || 2 where r est i (respectively, r gt i ) is the ray direction of the i-th pixel by the estimated (respectively, GT) multi-camera calibration. Our distance is d = √ e(R)/N where N is the number of (sampled) rays in a multi-camera image. Note that d is expressed in radians if d 1; we always convert it in pixels by dividing it by the angular resolution r in Tab. 1.

Frame-accurate synchronization

Here we experiment the FA synchronization summarized in Sec. 3 and detailed in Sec. 5. Fig. 7 draws the IAV (defined in Eq. 8) for consecutive frames taken in a rectilinear segment of the trajectory and the correlation function (that maps FA offset candidate o 0,1 to ZNCC 0,1 ) for cameras 0 and 1. This is done for biking (BC1), walking (WT) and car+mast (CC). There are similar variations of the IAV for different cameras and a single maximum of the ZNCC except for WT. Two consecutive offsets of WT have very similar greatest ZNCC values and the other ZNCC values are below, which suggest a half-frame residual time offset. These examples can convince the reader that we have enough information in the IAV to obtain a FA synchronization (at least if the cameras are helmet-held or mounted on a car thanks to a mast).

Tab. 2 shows the FA time offsets for all sequences. First we examine o j, j+1 for the four Gopro Hero 3 cameras in BC1, WT and FH. Different experiments have different o j, j+1 although they are taken by the same cameras. Thus synchronization should be done at every experiment. Furthermore, the wifi-based synchronization of the Gopro is not very accurate: about 0.04s and sometimes above 0.1s (reminder: their FpS is 100Hz or 48Hz). Second we check that the o j, j+1 are FA accurate when their ground truths are known (BC2, CC and WU).

Intrinsic parameters using GS.X.FA.INT

Before exploring BA for rolling shutter and synchronization in the next Sections, Sec. 8.4 experiments BA using two approximations: global shutter and zero sub-frame residual time offsets. In other words, we investigate the BA in [START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF] with two modifications: estimating the intrinsic parameters and minimizing the reprojection errors in the right image space. Sec. 8.4.1 compares error minimizations in the original and rectified image spaces using the polynomial distortion model. Sec. 8.4.2 compares results obtained without or with the central approximation.

We remind that the BA input is obtained as summarized in Sec. 3: initialization of multi-camera calibration on the video beginning (SfM and GS.C.FA.INT applied to the 2k first frames) followed by multi-camera SfM applied to the whole video. Note that the FoV angles used for equiangular initialization (Sec. 4) are in the FoV column of Tab. 1.

8.4.1.

Polynomial distortion model: original vs. rectified errors Tab. 3 compares the calibrations obtained by minimizing the reprojection errors in original (i.e. distorted) and rectified (i.e. undistorted) image spaces using GS.C.FA.INT applied to the 2k first frames of the CC video (the beginning of our method). Column "init" gives details on the initialization of the calibration: "72" means that we use an initial FoV angle equals to 360/5 • for monocular camera and "pat" means that we use the calibration estimated using a planar calibration pattern [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration ?[END_REF]. Although the numbers of 2D inliers are similar, the calibration error d of the distorted case is quite better (about 6 times smaller) than that of the rectified one. Such a difference can be explained as follows: there are large distortions between rectified and distorted images, the BAs minimize errors in different image spaces (rectified and distorted), and the distorted space is the right one to obtain a Maximum Likelihood Estimator (more details at the end of Sec. 2.4). Tab. 3 also provides the numbers of 2D inliers and RMS for GS.FA that enforces calibration "pat" during BA; our RMS and inliers are slightly better but the calibration error d of "pat" is the best. In the paper remainder, we always use the reprojection error in the original space. Second we compare GS.C.FA and GS.NC.FA and see that the non-central refinement (NC) changes almost nothing if INT is not refined. If INT is refined, the NC refinement improves the BC2 calibration but it does not improve (even degrades) the CC calibration. We interpret this result as follows: the central approximation is more tenable for CC than for BC2, since CC has a larger ratio between camera-scene distance and baseline than BC2. We also note that the RMS and 2D inliers are similar in all cases.

Tab. 5 provides accuracies of the intrinsic parameters of the first BC2 camera using GS.C.FA.INT and GS.NC.FA.INT. The absolute errors of f x , f y , u 0 , v 0 are about 2 pixels or less; the relative errors of k 1 is good and those of k i are bad if i > 2. The NC-values of f x , f y and u 0 are slightly better than those of C.

Rolling shutter and sub-frame accurate synchronization

First Sec. 8.5.1 provides all estimation errors of several BAs for videos BC2 and CC, that have complete ground truth. Second Sec. 8.5.2 provides SFA time offsets, line delay, and top view of reconstruction for every video using RS.C.SFA.INT.

Accuracies

Tab. 6 provides the errors e(∆), e(τ) and d (Sec. 8.2) for several BAs estimating both SFA time offsets ∆ j and line delay τ (Sec. 6). We compare GS. First we experiment on the only RS sequence that has complete ground truth: BC2. We see that the simultaneous estimation of INT and RS.SFA has quite smaller e(τ) and smaller d than separate estimations (both C and NC BAs). However the separate case has a twice smaller e(∆) than the simultaneous case, which in turn is more than twice smaller than that of the SFA refinement without BA. We remind that e(∆) cumulates absolute errors of SFA synchronization: a value of 0.1 (for the simultaneous case) means that the mean SFA sync. error of n cameras is only 0.1/(n -1). The NC BAs also greatly reduce d.

Second we experiment on CC, which is the only real sequence with complete ground truth. Since it is GS, the relative error e(τ) is not a number and we replace it by y max f τ (the smaller absolute value, the best result). We see that all BAs provides small |y max f τ| compared to that of consumer RS cameras: we obtain values in [0.002, 0.009], which are small compared to typical values in [0.8, 0.9] of consumer RS cameras. Furthermore, all e(∆) are smaller than 0.055 for five cameras; d increases and e(∆) decreases by the NC BAs. The SFA refinement without BA provides the smallest e(∆).

Time offsets, line delays and reconstruction

Tab. 7 shows time offsets f ∆ j , normalized line delay y max f τ and error e(τ) for all videos by applying RS.C.SFA.INT. We see that error e(τ) is less than 7.2% except in the WT case. In the WT case, τ is over-estimated (y max f τ is even greater than its theoretical maximum value 1) and has large error e(τ) equal to 16%. In contrast to this, e(τ) in the BC1 case looks lucky. In fact, the τ value in a single experiment should be moderated since it depends on the keyframe choice (this will be experimented in Sec. 8.6). Note that a negative τ (for WU) simply means that the time of the y-th line increases when y decreases.

Last Figs. drift is less noticeable in the other examples). There are several reasons: we do not enforce loop closure, the incremental multicamera SfM by local BA [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] is done using an intermediate calibration computed from only 2k first frames, and the final BA (RS.C.SFA.INT) does not remove the drift. We redo the incremental SfM using the final multi-camera calibration (computed from the whole sequence by RS.C.SFA.INT) and see that an important part of drift is removed. This suggests that the final multi-camera calibration is better than the intermediate one.

Stability with respect to keyframe sampling

Now we experiment the stability of our results (SFA synchronization, line delay and calibration) with respect to moderated changes of the keyframe sampling. The keyframe sampling is tuned by a single threshold N 3 , which is a lower bound for the number of matches between three consecutive keyframes (more details in Appendix G). For every value N 3 ∈ {400, 425, 450, 475, 500}, we apply multi-camera SfM based on keyframe sampling followed by RS.C.SFA.INT and then discuss the results. The initial multi-camera calibration and FA synchronization are the same for all N 3 and are computed from the video beginning as in the other experiments.

The left of Tab. 8 shows estimation errors e(∆), e(τ) and d. There are 206 keyframes if N 3 = 400 and 248 keyframes if N 3 = 500. The variations of errors are important: from single to double for e(∆), from single to quadruple for e(τ), and about 30% for d. We provide an explanation in Sec. 8.6.1 and a correction of the results in Sec. 8.6.2. Tab. 9 shows time offsets f ∆ j and error e(τ) for the longest sequence BC1. The variation of e(τ) are also important; the variation of f ∆ j are less than 0.032.

Analysis

We remind that the reprojection errors in the i-th keyframe are computed using the approximation in Eq. 12 of the multicamera trajectory M(t) where t ≈ t i : M(t) is a linear combination of m i-1 , m i and m i+1 , and M is linear in time tt i . The better this approximation, the smaller the reprojection errors in the i-th keyframe.

At first glance, the estimation errors decrease if N 3 increases: if N 3 increases, the keyframe density increases, thus the accuracy of these approximations is better (the remainders O(∆) and O(∆ 2 ) in Sec. 6.4 decrease), the reprojection errors decrease and last the estimation errors decreases). However, the errors in Tab. 8 are not decreasing series but look noisy.

Here is a second explanation. The true value of M(t) near t i can be different to its approximated value for some i, e.g. if the true speed vector M (t i ) is different to the speed vector D i 1 computed using the multi-camera poses of neighboring keyframes i -1 and i + 1 in Eq. 11. Then a keyframe with bad approximation has high reprojection errors and act as an outlier perturbing the BA. The estimation errors in the left of Tab. 8 depend on the set of this kind of outliers, which in turn depends on N 3 .

Correction

The idea is simple: if the i-th keyframe has high reprojection errors, we redefine its approximation by M(t) = m i + (tt i )d i if t ≈ t i thanks to a new velocity parameter d i ∈ R 6 that is estimated by BA like m i . Then the camera motion is not constrained by keyframes i -1 and i + 1 if t ≈ t i , and we expect that , all errors are improved in the following sense: both mean and maximum of every error are reduced, the variations of e(∆) and d are damped (the variations of e(τ) expressed using ratio max/min are not damped due to a small error 0.4% for N 3 = 500). Last Tab. 10 shows time offsets f ∆ j , normalized line delay y max f τ and error e(τ) for all videos by applying RS.C.SFA.INT+RS.C.SFA.INT.h. We see that all errors e(τ) have the same magnitude order (in interval [2.8, 8]). This contrasts to Tab. 7 using N 3 = 450, where e(τ) is small for BC1 and large for WT.

Name #d i f ∆ 1 f ∆ 2 f ∆ 3 y max f τ e(τ

Robustness with respect to FA synchronization

We would like to know whether an error in the FA synchronization can be corrected by SFA synchronization. Such an error can have two reasons: the IAV variations are insufficient in the video beginning for FA accurate synchronization, or a consumer camera skips frame(s) for any technical reasons. We simulate such an error by skipping x frames of camera 1 (reminder: camera 0 is the first one), then we use multi-camera SfM followed by RS.C.SFA.INT and compare the results for x ∈ {0, 1, 2, 3}.

Tab. 11 shows errors e(∆), e(τ), d and time offsets f ∆ j estimated for BC2. We see that e(∆) and d increases moderately if x = 1, e(∆) is multiplied by 2.5 and e(τ) by 3.9 if x = 2, and all errors increase a lot if x = 3.

Variations of τ and ∆ j in a long sequence

Here we examine the variations of τ and ∆ j in a long sequence. We split the GS.C.FA.INT reconstruction of BC1 into six segments of 300 keyframes (segments 0-299, 300-599, etc) and independently apply RS.C.SFA.INT to every segment. Tab. 12 shows the results. The variations of f ∆ j are moderated (less than 0.2) and the f ∆ j globally increase over time, i.e. when s increases. At first glance, we could expect that we can detect a frame skipped by a camera (if any) by an increase/decrease of 1 as in Sec. 8.7. However this is not the case. Since all cameras are in the same manufacturing series and have the same setting, all cameras skip similar numbers of frames (if any) in a segment, which in turn would imply that we do not observe large time offset perturbations. Furthermore we could interpret the slow increase of the f ∆ j s as follows: the FpS of camera 0 is slightly lower than those of the other cameras.

x f ∆ 1 f ∆ 2 f ∆ 3 e (
The variations of τ are important (especially in the first half part of BC1) and τ globally decreases over time. The variations of τ are reduced if we take a larger video segment size, e.g. we divide by two the y max f τ range (max-min in Tab. 12) with 500 keyframes per segment.

Other experiments

All previous experiments are done using the linear approximation of M(t) in Eq. 12 and using the approximation t p = t p in Eq. 17. Tab. 13 shows the results for all videos using the quadratic approximation of M(t) in Eq. 15 or using Sec. 7.3 (i.e. without the approximation t p = t p in Eq. 17). We do not observe significant improvements of e(τ) compared to those in Tab. 7 except for FH using Eq. 15. By recomputing errors e(∆) and d for BC2 and CC using these two changes, we obtain very similar results as in Tab. 6: e(∆) difference is less than 0.004 and d difference is less than 0.01.

Conclusion

This article introduces the first self-calibration method for a multi-camera moving in an scene, that simultaneously estimates intrinsic parameters, inter-camera poses, time offsets and line delay in addition to the usual parameters (3D points and multicamera poses). We start by a rough calibration assuming that the multi-camera is central and omnidirectional without privileged direction. Then we estimate frame-accurate time offsets using monocular structure-from-motion and bundle adjustment (SfM and BA) without assumption on the field-of-view shared by adjacent cameras. Last we apply multi-camera SfM and BA twice: using simple and complicated camera models. The former forces to 0 line delay, sub-frame residual time offsets and baseline between cameras; then the former initializes the latter.

We experiment in a context that we believe useful for applications (3D modeling and 360 videos): several consumer cameras or a spherical camera mounted on a helmet and moving along long trajectories by walking and biking (among others). Long trajectories are useful for calibration accuracy and are allowed since our BA only refines the keyframes provided by SfM. We compare central and non-central results, provide accuracy for calibration/line delay/time offsets with respect to ground truth, examine the influence of the tuning of keyframe selection, show variations of time offsets in a long sequence, experiment and compare different approximations (for time continuous camera trajectory and for reprojection errors).

However the method has a limitation: the image deformations due to rolling shutter should be moderated for SfM (before the final BA that refines all parameters including line delay) since SfM uses the global shutter approximation. Furthermore several improvements and future work are possible. First the initialization, which is not the paper topic, can be improved thanks to previous work for both intrinsic parameters and inter-camera poses. Second a preprocessing should select segment(s) in the video where we safely apply SfM. Third variants of the method can be experimented, e.g. by using nonminimal parametrization of rotations, alternative keyframe selections and other camera models. Last, we should examine the improvements in applications provided by our non-zero line delay and sub-frame-accurate time offsets. Second we use Eq. C.1 with t = t i , a = t i+1t i , b = t it i-1 , ∆ = max i (t i+1t i ), and obtain M (t i ) = D 1 (M(t i-1 ), M(t i ), M(t i+1 ), t i+1t i , t it i-1 ) + O(∆ 2 ). In practice, it is improbable that a product or sum of nonzero matrices is zero. Thus we replace the inclusions above by equalities in our proof, i.e. we use L(A B) = L(A) × L(B) and L(C) = ∪ i L(C i ).

(E.8)

Let Z g and Z r be the 6m × 6m top-left blocks of the RCS in the standard case and our case (more details in Sec. 6.7.2), respectively. In the next section, we show that L 6×6 (Z r ) = L 6×6 (Z g ) + {-1, 0, 1} 2 .

(E.9)

In these expressions and the following ones, we implicitly omit integers that are below 1 and above m (e.g. we omit i -1 if i = 1 and omit i + 1 if i = m).

Appendix E.2. Proof of Eq. E.9

The image projection function of the l-th 3D point in the ith multi-camera pose is ϕ g i,l in the standard case and ϕ r i,l in our case. According to Sec. 6.7.2, we have ∂x l and WV -1 W = l W l V -1 l W l , we have L 6×6 (WV -1 (W) ) = ∪ l (L 6×3 (W l )) 2 . Thus L 6×6 (W g V -1 (W g ) ) = ∪ l (V l ) 2 and (E.20)

L 6×6 (W r V -1 (W r ) ) = ∪ l (V l + {-1, 0, +1}) 2 (E.21) = {-1, 0, +1} 2 + ∪ l (V l ) 2 .
(E.22) Thus L 6×6 (Z g ) = ∪ i {(i, i)} ∪ l (V l ) 2 and (E.23)

L 6×6 (Z r ) = {-1, 0, +1} 2 + (∪ i {(i, i)} ∪ l (V l ) 2 ). (E.24)
We obtain Eq. E.9.

Appendix F. Proof of Eq. 29 Last we obtain Eq. 29 using Eqs. F.7, F.8 and 23.

Figure 1 :

 1 Figure 1: Time continuous trajectory of a multi-camera. Left: four monocular cameras at time t i , which have non-zero time offsets. Right: a rolling shutter monocular camera, which moves and sees points at several times/lines in a single frame.

Figure 2 :

 2 Figure 2: Shape of the standard (left) and our (right) Z for a video sequence with 1573 keyframes and closed loops and same inliers.

Figure 3 :

 3 Figure 3: Cameras (four Gopro Hero 3 in a cardboard) and images for BC1, WT and BC2. The rolling shutter always goes from right to left, the image motion goes toward left on the two left columns and goes toward right on the two right columns.

  FpS f , angular resolution r (millirad.), diameter b of multi-camera centers, line delay τ (ground truth), time offset f ∆ j (ground truth), approximate trajectory length l, numbers of frames f r and keyframes k f r, FoV angle used for initialization, maximum of angles |β i | (radians) of our parametrization in Sec. 6.6.2.

Figure 4 :

 4 Figure 4: Cameras (four Gopro Hero 3 in their housings) and images for FH. The rolling shutter always goes from top to bottom, the image motion goes forward/toward right/backward/toward left.

Figure 5 :

 5 Figure 5: Spherical camera (Ricoh Theta S) and image for WU. Both rolling shutter and image motion go from bottom to top (τ is negative, the bottom line is the oldest).

Figure 6 :

 6 Figure 6: Cameras (PointGrey Ladybug 2) and images for CC. The cameras are global shutter. The camera pointing toward the sky is not used in experiments.

Figure 7 :

 7 Figure7: IAV for two cameras (left) and their correlation curves (right) for rectilinear trajectory segments of trajectories of BC1 (top), WT (middle) and CC (bottom). In the left, we have frame numbers (x-axes) and IAVs in radian (y-axes). In the right, we have offset candidates o 0,1 ∈ Z (x-axes) and its correlation ZNCC 0,1 ∈ [-1, 1] (y-axes).

8.4. 2 .

 2 Central vs. non-central Tab. 4 compares calibration error d (and RMS and 2D inliers) obtained using GS.C.FA.INT, GS.NC.FA.INT, GS.C.FA and GS.NC.FA applied to videos BC2 and CC. First GS.C.FA.INT provides a better (smaller) d than GS.C.FA since the intrinsic parameters (INT) are estimated from a longer video. The comparison between GS.NC.FA.INT and GS.NC.FA is similar.

Figure 8 :

 8 Figure 8: Top views of RS.C.SFA.INT reconstructions of BC1, WT, FH, BC2 and CC (from top to bottom) without loop closure. The input videos of BC1, WT and FH are taken by four Gopro cameras mounted on a helmet. The FH trajectory has a lot of sharp S turns.

Figure 9 :

 9 Figure 9: Top views of WU reconstructions without loop closure. The input video is taken by the Ricoh Theta S mounted on a helmet. The drift is between the two arrows. Top: result of RS.C.SFA.INT. Bottom: incremental SfM [28] that is redone using the calibration estimated by RS.C.SFA.INT.

Figure A. 10 :

 10 Figure A.10: Notations in two cases: ξ > 1 (left) and 0 < ξ < 1 (right).

(C. 6 )

 6 If a matrix A is partitioned by blocks A i, j (both horizontally and vertically), we define a list of index pairsL(A) = {(i, j), A i, j 0}. (E.5)If matrices A and B are horizontally partitioned by blocks A i and B j respectively, A B is partitioned by blocks A i B j and for these blocks we haveL(A B) ⊆ L(A) × L(B). (E.6) If C = i C i and all C i have the same block partition, L(C) ⊆ ∪ i L(C i ). (E.7)If all blocks considered by L have the same size a × b, we can write L a×b instead of L.

  ϕ g i,l = ϕ(m i , m , x l ) and ϕ r i,l = ϕ(m i-1 , m i , m i+1 , m , x l ). (E.10) i ∈ {i -1, i, i + 1}. We rewrite this using the notations in Appendix E.1:) = {i -1, i, i + 1}. (E.11)Notations ϕ i,l and Z, U • • • are used in expressions that hold for both "r" and "g" upper-indices added to these notations. Let V l be the list of keyframe indices where the l-th point is inlier. According to Eqs. 21 and 22, we have Z = U -WV -1 W where

Letθ=,

  ( f x , f y , u 0 , v 0 , k 1 , k 2 , • • • , k n , x j ), zu = π(x j ), (F.2) g(z, θ) = (1 + n i=1 k i ||z|| 2i )z -zu . (F.3)Thanks to Eqs. 1 and 2 and since p = z d , we have g(p, θ) = 0. First we note that∂z ∂(u 0 , v 0 ) f x , f y ) = ∂g ∂z ∂z ∂( f x , f y )

Table 1 :

 1 Datasets:

			3 100	1.56	7.5	9.10	?	2500 50.4k	90	2047	343k	0.223
	WalkTown (WT)	4*Gopro 3 100	1.56	7.5	9.10	?	900 70.3k	90	1363	240k	0.268
	FlyHill (FH)	4*Gopro 3 48	1.06	18	11.3	?	1250 8.6k	90	627	432k	0.494
	BikeCity2 (BC2)	4*Gopro 3 100	1.56	7.5	9.12	i/4	615 12.5k	90	225	51k	0.074
	CarCity (CC)	Ladybug 2 15	1.90	6	0	0	2500 7.7k	72	891	282k	0.068
	WalkUniv (WU)	Theta S	30	3.85	1.5	-32.1	0	1260 29.4k 200	1287	154k	0.129

Table 2 :

 2 o 0,1 o 1,2 o 2,3 o 3,4 Zncc 1 Zncc 2 FA time offsets o j, j+1 with loop constraint. Zncc 1 is the greatest sum of the ZNCCs of the n computed time offsets, and Zncc 2 is the second greatest ZNCC (thus -n ≤ Zncc i ≤ +n). We remind that o j, j+1 counts a signed number of frames between the j-th and the j + 1-th videos (it is "na" if the multi-camera has less than j + 2 cameras, i.e. if n < j + 2).

	BC1 -5	3	4	na	3.912 3.884
	WT -15	-1	14	na	3.919 3.918
	FH	-1	1	-2	na	3.991 3.983
	BC2	0	0	0	na	3.915 3.907
	CC	0	0	0	0	4.987 4.347
	WU	0	na	na	na	0.993 0.677

Table 3 :

 3 Comparing accuracy of gs.fa.X using rectified and distorted reprojection errors on 2k first frames of CC (reminder: d is in pixels).

	Error	init	Method	#2D		RMS	d
	Rectified 72 gs.c.fa.int 213335 1.216 9.575
		pat	gs.fa	213015 1.225 1.023
	Distorted 72 gs.c.fa.int 213495 0.932 1.683
		pat	gs.fa	213108 0.946 1.023
	Method	d	RMS #2D	d	RMS #2D
	gs.c.fa	3.379 0.728 204k 1.685 0.938 965k
	gs.c.fa.int	2.018 0.723 204k 1.173 0.938 965k
	gs.nc.fa	3.397 0.727 204k 1.684 0.938 965k
	gs.nc.fa.int 1.417 0.723 204k 1.313 0.937 965k

Table 4 :

 4 Accuracies of gs.fa.X for BC2 (left) and CC (right).

Table 5 :

 5 C.FA.INT+RS.C.SFA and Accuracies of intrinsic parameters of the first camera (BC2) using gs.X.fa.int. Reminder: we use absolute errors in pixels for f x , f y , u 0 , v 0 and relative errors for the k i s.RS.C.SFA.INT, i.e. we compare separate and simultaneous estimations of INT and RS.SFA parameters. We also compare these central BAs and their non-central versions, and the SFA synchronization without BA in Sec. 5.2.

			gs.c.fa.int	gs.nc.fa.int
		G.T.	value	error	value	error
	f x 580.773 582.809 2.037 581.296 0.523
	f y	581.266 582.844 1.578 582.605 1.339
	u 0 640.827 640.989 0.162 640.978 0.151
	v 0 469.056 471.305 2.249 471.540 2.284
	k 1	0.368	0.368	7e-4	0.368	6e-4
	k 2	0.067	0.063	0.061	0.062	0.062
	k 3	0.013	0.026	1.001	0.025	0.970
	k 4	0.002	-0.013	6.463	-0.013	6.378
	k 5	0.013	0.018	0.434	0.018	0.420

Table 6 :

 6 8 and 9 show a top view of the RS.C.SFA.INT reconstructions (both keyframes locations and 3D point cloud).In the WU case, we observe a non-negligible drift since the beginning and end of the trajectory should be the same (the Accuracies of rs.sfa.X.(int) for BC2 and for CC.

	Methods applied to BC2 e(∆)		e(τ)	d
		gs.c.fa.int+rs.c.sfa	0.057 14.6% 1.970
		rs.c.sfa.int	0.097	2.7%	1.476
		gs.nc.fa.int+rs.nc.sfa	0.051 12.2% 1.312
		rs.nc.sfa.int	0.111	3.7%	0.366
		gs.sfa (in Sec. 5.2)	0.215	na	na
		Methods applied to CC	e(∆)		y max f τ	d
		gs.c.fa.int+rs.c.sfa	0.052 -0.0052 1.176
		rs.c.sfa.int	0.055 -0.0069 1.167
		gs.nc.fa.int+rs.nc.sfa	0.034 -0.0020 1.322
		rs.nc.sfa.int	0.039 0.0086 1.313
		gs.sfa (in Sec. 5.2)		6e-3		na	na
		f ∆ 1	f ∆ 2	f ∆ 3		y max f τ	GT	e(τ)
	BC1 -0.334 -0.153 0.132	0.8755	0.8736	0.2%
	WT -0.583 -0.320 -0.795		1.013	0.8736 16.0%
	FH	0.287	0.203 -0.326 0.8372	0.7810	7.2%
	BC2 0.246	0.546	0.797	0.8989	0.8755	2.7%
	CC -0.017 -0.013 -0.006 -0.0069	0	nan
	WU 0.001	na	na	-0.8882 -0.9244 3.9%

Table 7 :

 7 SFA time offsets and line delay accuracy for all datasets using rs.c.sfa.int. Here GT is the ground truth of y max f τ.

Table 8 :

 8 Accuracy stability with respect to keyframe sampling for BC2 using rs.c.sfa.int (left) and rs.c.sfa.int+rs.c.sfa.int.h (right) with h = 70%.

	N 3	k f r	f ∆ 1	f ∆ 2	f ∆ 3	y max f τ e(τ)
	400	1813 -0.341 -0.171 0.131 0.8920 2.1%
	425	1929 -0.337 -0.155 0.131 0.8882 1.6%
	450	2047 -0.334 -0.153 0.132 0.8755 0.2%
	475	2166 -0.366 -0.156 0.134 0.9399 7.5%
	500	2256 -0.337 -0.141 0.130 0.8786 0.5%
	mean	2042 -0.343 -0.155 0.132 0.8948 2.4%
	mx-mn 443	0.032	0.030 0.004 0.0644 7.3%

Table 9 :

 9 Stabilities of time offsets and line delay with respect to keyframe sampling for BC1 using rs.c.sfa.int. The number of keyframe is k f r.

Table 10 :

 10 SFA time offsets and line delay accuracy for all datasets using rs.c.sfa.int+rs.c.sfa.int.h with h = 70% (to be compared with Tab. 7). The number of keyframes with additional velocity parameter is #d i .the resulting reprojection errors of the i-th keyframe decrease such that the i-th keyframe does no act as an outlier of the BA.In practice, we start from a current estimation obtained by RS.C.SFA.INT and introduce an user defined percentage h. Let c h be the h-fractile over all reprojection errors. For every keyframe, we compute the RMS of its own reprojection errors. If this RMS is greater than c h , the keyframe has an additional velocity parameter d i as above (otherwise it does not have). We name RS.C.SFA.INT.h the new BA obtained by modifying RS.C.SFA.INT like this. Tab. 8 provides estimation errors e(∆), e(τ), d of both RS.C.SFA.INT and RS.C.SFA.INT+RS.C.SFA.INT.h using h = 70%. Thanks to the correction

Table 11 :

 11 Accuracies of rs.c.sfa.int applied to BC2 if we skip x additional frame(s) of camera 1. The ideal result meets f ∆ 1 = x + 0.25.

				∆)	e(τ)	d
	0 0.246 0.546 0.797 0.097 2.7% 1.476
	1 1.342 0.505 0.769 0.117 1.0% 1.743
	2 2.026 0.489 0.760 0.245 10.5% 1.298
	3 2.472 0.532 0.781 0.841 10.4% 3.244
	s	f 1	f ∆ 2	f ∆ 3	y max f τ	e(τ)
	0	-0.443 -0.180 0.014 1.0570 20.1%
	1	-0.417 -0.152 0.014 0.9347 6.9%
	2	-0.306 -0.073 0.028 0.7678 12.2%
	3	-0.308 -0.163 0.115 0.9225 5.6%
	4	-0.271 -0.157 0.178 0.9056 3.7%
	5	-0.299 -0.105 0.199 0.8394 3.9%
	mean	-0.341 -0.138 0.091 0.9045 8.7%
	max-min 0.172	0.107 0.185 0.2892 16.4%

Table 12 :

 12 Stabilities of time offsets and line delay over time in long sequence BC1 using rs.c.sfa.int. The s-th video segment is taken between keyframes 300s and 300(s + 1) -1.

Table 13 :

 13 SFA time offsets and line delay accuracy for all datasets using rs.c.sfa.int (to be compared with Tab. 7).

  6×6 (U) = ∪ i∈V l (L 2×6 ( ∂ϕ i,l ∂M )) 2 , we have L 6×6 (U g ) = ∪ i {(i, i)} and (E.14) L 6×6 (U r ) = ∪ i {i -1, i, i + 1} 2 (E.15) = {-1, 0, +1} 2 + ∪ i {(i, i)}. (E.16)Furthermore, W in Eq. 21 is horizontally partitioned in blocks W l = i∈V l (

	U = Since L ∂ϕ i,l i∈V l ( ∂ϕ i,l ∂M ) ∂ϕ i,l ∂M , W = i∈V l ( ∂ϕ i,l ∂M ) ∂ϕ i,l ∂X , (E.12) V = i∈V l ( ∂ϕ i,l ∂X ) ∂ϕ i,l ∂X . (E.13) ∂M ) L 6×3 (W l ) = ∪ i∈V l L 6×3 (( ∂ϕ i,l ∂x l . Since = ∪ i∈V l L 2×6 ( ∂ϕ i,l ∂ϕ i,l ∂M ∂M ), ) we have L 6×3 (W g l ) = V l and L 6×3 (W r l ) = V l + {-1, 0, +1}. ∂ϕ i,l ∂x l ) Since V in Eq. 21 is block-wise diagonal with invertible blocks (E.17) (E.18) (E.19) V l = i∈V l ( ∂ϕ i,l ∂x l ) ∂ϕ i,l

Preprint submitted to Computer Vision and Image Understanding December 7, 2017

Acknowledgments

Thanks to CNRS, Université Clermont Auvergne and Institut Pascal for funding Maxime Lhuillier and Thanh-Tin Nguyen.

Appendix A. Properties of the unified camera model

We remind that this model is described in Sec. 4.2 using a perspective camera. Let k = 0 0 1 , c = 0 0 -ξ and 3D point x. Let ν be the angle between k and a ray of the perspective camera (a half-line started at c and including x/||x||). Let µ be the angle between k and a ray of the unified camera (a half-line started at 0 with direction x/||x||).

Appendix A.1. Theoretical field-of-view Fig. A.10 shows notations ν, µ, µ 0 (a value of µ), x/||x||, c and k in two cases: ξ > 1 and 0 < ξ < 1. In both cases, µ 0 is the maximum value of µ that ensures that there is only one backprojected ray direction corresponding to the projection p(x) by the unified camera model. The angle µ 0 is the half-angle of the theoretical FoV of this model, i.e. by ignoring the bounded size of the image and the projection of the camera itself. The FoV in the paper core ignores nothing and is included in the theoretical FoV.

If ξ = 0, the unified camera model is the standard perspective camera model with the camera center c = 0. Since x/||x|| is in front of the perspective camera, the theoretical FoV is the halfspace z > 0 and µ 0 = π/2.

If 0 < ξ < 1, c is inside the unit sphere S. Since x/||x|| is in front of the perspective camera, the theoretical FoV is the half-space z > -ξ and cos µ 0 = -ξ.

If ξ > 1, c is outside S and S is entirely in front of the perspective camera. The projection of S by p is an ellipse and its interior. Let C be the cone that is tangent to S with apex c, i.e. the union of every line (cy) that intersects S at a single point y. We have y (c-y) = 0 and k y < 0 and cos µ 0 = k y = 1/(-ξ). For example, ξ = 2 (in Sec. 4.2.2) implies that µ ≤ µ 0 = 2π/3.

Appendix A.2. Angle of back-projected ray

There is a relation between ν and µ in all cases. Using notation x 1 x 2 x 3 = x/||x||, we have cos µ = x 3 and

Since p(x) = ( f x x 1 /(ξ + x 3 ) + u 0 , f y x 2 /(ξ + x 3 ) + v 0 ) and p(k) = (u 0 , v 0 ) = z 0 ,

Appendix B. IAV property

Here we show that the IAV in Eq. 8 is the same for two frames of different but jointly moving cameras if they are taken at the same time. We remind properties of a change of basis in R 3 expressed by rotation matrices: the columns of R A,B are vectors of B expressed using coordinates in A, we have

Since the monocular SfMs are not done in the same coordinate system, the notation R t i used in Eq. 8 is ambiguous. Here we write instead R t w i ,i where w i is the (world) basis where is reconstructed the i-th video (vectors of the i-th camera at frame t expressed using coordinates in w i ). Furthermore, we note that R w i ,w j and R i, j do not depend on frame numbers (the former is obvious, the latter is due to the fact that the cameras are rigidly mounted). If the t i -th frame of the i-th camera and the t j -th frame of the j-th camera are taken at the same time,

Since the cameras have the same FpS, we also have

Thanks to Eqs. B.1 and B.2, we obtain

w j , j (R

Since trace(XY) = trace(YX), we obtain trace(R

Thus θ t j j = θ t i i according to Eq. 8.

Appendix C. Derivatives of the camera motion M(t)

The notations used in Appendix C.1 and Appendix C.2 are defined in Secs. 6.4.1 and 6.4.2, respectively.

First we show that

First we show that

We eliminate M (t) by summing b(Eq. C.2)+a(Eq. C.3):

Since a > 0 and b > 0,

Second we use Eq. C.7 with t = t i , a = t i+1t i , b = t it i-1 , ∆ = max i (t i+1t i ), and obtain

Appendix D. Estimation of rotations A and B

Here we not only compute A and B in the definition of our rotation parametrization R (Eq. 20), but also initialize E M (t i ). Let R 0 i be the rotation of the i-th keyframe estimated by GS BA (before the final BA in Sec. 6). Using the definitions of E M and t i in Secs. 6.2 and 6.3,

Since the multi-camera trajectory has small pitch and roll, there are rotations A and B such that ∀i, A -1 R 0 i B -1 is almost a rotation around k (details in Appendix D.1). Let angles (α i , β i , γ i ) be such that E(α i , β i , γ i ) = A -1 R 0 i B -1 and E is defined in Eq. 19 and β i is close to 0 (details in Appendix D.2). We initialize E M (t i ) = α i β i γ i . Thanks to Eqs. 20 and 19, we obtain

According to Sec. 6.6.2, E M (t i ) is a singularity of E iff β i ∈ π/2 + πZ. Since R has the same singularities as E (proof in Appendix D.3) and β i is close to 0, E M (t i ) is far from the singularities of R.

Appendix D.1. Technical Details: Estimate A and B Let R(v, θ) be the rotation with axis v and angle θ. Since the multi-camera trajectory has small pitch and roll (Sec. 6.6.2), all (R 0 i ) T R 0 j are roughly rotations sharing a same axis v ∈ R 3 . Thus there are rotation R and angle θ i such that R 0 i ≈ RR(v, θ i ). For all i and j, (R 0 i ) T R 0 j ≈ R(v, θ j -θ i ). Let v i, j be the axis of (R 0 i ) T R 0 j . First we search v as the most colinear vector to all v i, j , i.e. v maximizes i, j (v T i, j v) 2 . Thus v is the eigen vector of the largest eigen value of the symmetric matrix i, j v i, j v T i, j . Second we estimate rotation R such that RR 0

Since E is surjective on the set of 3D rotations, the angles α i , β i and γ i exist. Furthermore, they are defined up 2π multiples. We choose β i that has the smallest |β i |. Since E(α i , β i , γ i ) ≈ R(k, γ i ), β i is close to 0. We also remind that E M (t) is continuous (Sec. 6.2) and |t it i+1 | is small thanks to the keyframe sampling. Thus the γ i series is chosen such that |γ i -γ i-1 | is as small as possible, and we do similarly for α i (|β i -β i-1 | is also small).

Appendix D.3. R and E have the same singularities

Here we show that ker ∂R = ker ∂(ARB) if A and B are two invertible 3 × 3 matrices. Let x ∈ ker ∂R, R i j be the coefficients of R, and ∂R i j be the gradient of R i j with respect to parameters (α, β, γ). Thus ∂R i j .x = 0 and (∂(ARB) i j ).x = (∂(

We see that ∂(ARB).x = 0, i.e. ker ∂R ⊆ ker ∂(ARB). Appendix E.1. Notations and prerequisites If L 1 and L 2 are two lists of integers, we use notations

We also use Eq. E.1 if L 1 and L 2 are two lists of integer pairs. If a matrix A is partitioned by blocks A i (horizontally or vertically), we define an index list

Keyframe sub-sampling of the videos

The keyframe sub-sampling is that in Sec. 2.3 of [START_REF] Royer | Monocular vision for mobile robot localization and autonomous navigation[END_REF] with an improvement if the multi-camera motion is slow. It is based on counting the number of Harris points matched using ZNCC correlation between the current frame and its two preceding keyframes. Every frame is considered in increasing order of its index. First the current frame is rejected if its image motion compared to that of the last keyframe is small, i.e. if 70% of its matches have a 2D motion less than 5 pixels. Second, a non-rejected current frame is selected as a keyframe if ( 1) it has at least N 2 matches with the previous keyframe and ( 2) it has at least N 3 matches with the two previous keyframes and (3) the next frame does not meet (1) or (2). We always use N 3 = N 2 /2, use N 3 = 450 for the multi-frame of four Gopro cameras at 100 FpS, N 3 = 1200 for four Gopro cameras at 48 FpS, N 3 = 625 for Ladybug 2 and N 3 = 300 for Theta S.