Thanh-Tin Nguyen

Maxime Lhuillier

Adding Synchronization and Rolling Shutter in Multi-Camera Bundle Adjustment

Multi-cameras built by fixing together several consumer cameras become popular and are convenient for applications like 360 videos. However, their self-calibration is not easy since they are composed of several unsynchronized and rolling shutter cameras. This paper introduces a new bundle adjustment for these multi-cameras that estimates not only the usual parameters (camera poses and 3D points) but also the synchronization and the rolling shutter of the cameras. We experiment using videos taken by GoPro cameras mounted on a helmet, moving along trajectories of several hundreds of meters or kilometers, and compare our results to ground truth.

Introduction

Multi-cameras built by fixing together several consumer cameras become popular thanks to their prices, high resolutions, their growing applications including 360 videos (e.g. in YouTube) and generation of virtual reality content [1,2]. However such a multi-camera also has drawbacks. A first problem is the lack of accurate synchronization between the cameras. In the usual cases (e.g. GoPro), the camera manufacturer provides a wifi-based synchronization: the user starts all videos by a single click. However it is too inaccurate for applications such as 360 video and 3D modeling. Secondly, a low price camera implies that the camera is rolling shutter or RS. This means that two different lines of pixels of a frame are acquired at different instants (in a global shutter or GS camera, all pixels of a frame have the same time). Both inaccurate synchronization and RS complicate the self-calibration for the same reason: they act as time varying relative pose between the cameras, i.e. the multicamera has a varying non-central calibration [START_REF] Geyer | Geometric models of rolling-shutter cameras[END_REF].

This paper introduces a new bundle adjustment (or BA) for these multi-cameras, that simultaneously estimates not only the usual parameters (camera poses and 3D points) but also the synchronization and the RS coefficient. We start from an initial calibration with the simplest camera model (GS) and a frame-accurate synchronization (FA) provided by previous self-calibration methods [START_REF] Lhuillier | Synchronization and self-calibration for helmet-held consumer cameras, applications to immersive 3d modeling and 360 videos[END_REF]. FA means that we skip the first frames in the videos such that the sequels of the videos have the following property: the frames with the same frame index are taken at the same time up to the inverse of fps (frames-per-second). Our BA provides subframe-accurate synchronization (SFA), i.e. it estimates the residual time offsets between a reference camera and the others. It also estimate the RS coefficient, i.e. the time delay between two adjacent lines of a frame.

Previous work

In contrast to our multi-camera BA, the previous ones estimate neither synchronization nor RS. They [START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF][START_REF] Lhuillier | Synchronization and self-calibration for helmet-held consumer cameras, applications to immersive 3d modeling and 360 videos[END_REF][START_REF] Schneider | Bundle adjustment and system calibration with points at infinity for omnidirectional cameras[END_REF] assume that the cameras are synchronized and GS, only [START_REF] Klingner | Street view motion-from-structure-frommotion[END_REF] deals with known RS but needs other sensors. Previous monocular BA estimates 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] Duchamp | Multiple view 3D reconstruction with rolling shutter cameras[END_REF][START_REF] Hedborg | Rolling shutter bundle adjustment[END_REF]. In the context of visual SLAM [START_REF] Klein | Parallel tracking and mapping on a camera phone[END_REF], a GS BA is applied to a RS (monocular) camera thanks to a RS compensation: this method corrects beforehand the RS effects on the feature tracks by estimating instantaneous velocities of the camera.

Each 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. An assumption on the inter-frame motion is required. The BA in [START_REF] Duchamp | Multiple view 3D reconstruction with rolling shutter cameras[END_REF] adds extra parameters to avoid this assumption: it not only optimizes a pose but also rotational and translation speeds at every keyframe (not all frames). The translation speed is optional if its RS effect is negligible compared to that of rotation speed. 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-from-structure-frommotion[END_REF], the relative pose between an inter-frame pose and an optimized frame pose is provided by IMU at high frequency. The visual only RS approaches [START_REF] Duchamp | Multiple view 3D reconstruction with rolling shutter cameras[END_REF][START_REF] Hedborg | Rolling shutter bundle adjustment[END_REF][START_REF] Klein | Parallel tracking and mapping on a camera phone[END_REF][START_REF] Oth | Rolling shutter camera calibration[END_REF] are experimented on few meters long camera trajectories. Ours is experimented on longer trajectories (hundreds of meters, kilometers) since it only estimates poses at keyframes.

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 continuous-time 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 Structure-from-Motion (SfM). This is done thanks to an assumption on the multi-camera motion, which is tenable for helmet-held cameras in most cases.

3 Proposed Method

Initialization

First we assume that the monocular videos are approximately synchronized by removing few frames at their beginning (FA synchronization), i.e. the videos are synchronized up to the inverse of fps. The fps is assumed to be the same for all videos. 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 frames. Last we use standard SfM based on keyframe subsampling of the video and local BA [START_REF] Mouragnon | Generic and realtime structure from motion[END_REF], followed by global BA. The camera is self-calibrated assuming GS [START_REF] Lhuillier | Synchronization and self-calibration for helmet-held consumer cameras, applications to immersive 3d modeling and 360 videos[END_REF]. 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 function that maps Ω ⊆ R k to the set SO(3) of rotations of R 3 . We assume that there is a C 3 continuous function M : [0, 1] ⊂ R → R 3 × Ω that parametrizes the motion of the multi-camera. More precisely, M(t) T = T M (t) T E M (t) T where T M (t) ∈ R 3 is the translation and R(E M (t)) ∈ SO(3) is the rotation. The columns of matrix R(E M (t)) and T M (t) are the vectors of the multi-camera coordinate system at time t expressed in world coordinates. The choice of R is detailed in Sec. 4 for the paper clarity.

Thanks to these notations and assumptions, we will approximate M at every time t ∈ [0, 1] by using values of M at few times t 0 ,t 1 , • • •t n where t 0 = 0,t i < t i+1 and t n = 1. The M(t i ) are the only parameters of the multi-camera trajectory estimated by our BA. Sec. 3.3 defines t i and Sec. 3.4 describes our approximation of M(t) by using the M(t i ).

Time, Rolling-Shutter and Synchronization Parameters

The i-th keyframe is 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 fps. Since the cameras are RS, 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 time offset between the j-th camera and the 0-th camera: 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τ.

Approximations for the Multi-Camera Trajectory

First we have Taylor's expansion

M(t) = M(t i ) + (t -t i )M (t i ) + O(|t -t i | 2 ).
(1)

Second we provide a relation between derivative M (t i ) and all M(t i ). Let D be function

D(x, y, z, a, b) = bz a(a + b) - ax b(a + b) + (a -b)y ab , x, y, z ∈ R k+3 , a > 0, b > 0. ( 2
)
Let ∆ = max i (t i+1t i ) and shortened notation m i = M(t i ). Thanks to a linear combination of Taylor expansions of M at t i (more details in the supplementary material),

M (t i ) = D(m i-1 , m i , m i+1 ,t i+1 -t i ,t i -t i-1 ) + O(∆ 2 ). (3)
Third we approximate M(t) from the m i by neglecting all remainders expressed by "O" above. 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. 3.3). If 0 < i < n, M(t) = m i + (t -t i )D(m i-1 , m i , m i+1 ,t i+1 -t i ,t i -t i-1 ). ( 4
) If i = 0 (similarly if i = n), we use M(t) = m 0 + (t -t 0 )(m 1 -m 0 )/(t 1 -t 0 ).
Last we provide conditions that reduce the remainders of our Taylor developments, i.e. O(|tt i | 2 ) and O(∆ 2 ). If the density of keyframes in the video increases, ∆ decreases. Furthermore, |yτ| ≤ 1/ f ps and the FA synchronization provides |∆ j | ≤ 1/ f ps. Thus |t -t i | = |∆ j + yτ| decreases if the fps increases. If M = 0, both remainders are exactly 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 3D point x ∈ R 3 (in world coordinates) and its inlier observation p ∈ R 2 in the j-th sub-image 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 pp. Let (R j , t j ) be the pose of the j-th camera in the multi-camera frame. Let K j : R 3 \ {0} → R 2 be the projection function of the j-th camera. We assume that K 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 + ỹτ.

Second we detail the relation between p and x. Both E M (t p ) and T M (t p ), i.e. M(t p ), are defined by Eq. 4 using index i of the keyframe and t = t p . The coordinates of x in the multicamera coordinate system is x M = R(E M (t p )) (x -T M (t p )). The coordinates of x in the j-th camera coordinate system is x j = R j (x M -t j ). We also have p = K j (x j ).

Third we estimate p. We see that p needs the computation of x M , which in turn needs the computation of (the y coordinate of) p. This problem is solved thanks to an approximation in [START_REF] Klingner | Street view motion-from-structure-frommotion[END_REF]: t p is replaced by t p in the expression of x M , i.e. we assume that the multi-camera pose is the same at times t p and t p . We think that this is acceptable since |t p -t p | ≤ τ||p -p|| 2 and the magnitude order of τ is 10 -5 s/pixel and p is an inlier (i.e. ||p -p|| 2 ≤ 4 pixels).

Summary

First keyframes are selected and standard BA initializes the 3D points and the keyframe poses assuming GS cameras and FA synchronization (Sec. 3.1). Then the m i are computed from the poses. We also have τ = 0 (GS assumption) and ∆ j = 0 (FA assumption). Last we apply our BA which refines not only the m i and the 3D points, but also line delay τ (RS assumption) and/or the time offsets ∆ j (SFA assumption). It is based on Levenberg-Marquardt method.

Parametrization of Rotations

According to Sec. 3.2, R is a C 1 continuous function such that R(Ω) = SO(3) and Ω ⊆ R k . Following [START_REF] Oth | Rolling shutter camera calibration[END_REF], we prefer a minimal (non-redundant) parametrization R to avoid any constraints on the R entry and limit the number of estimated parameters. Sec. 4.1 is a reminder of these parametrizations and Sec. [START_REF] Furgale | Unified temporal and spatial calibration for multi-sensor systems[END_REF].2 details our choice.

Minimal Parametrizations of SO(3) for BA

During BA, the set of all rotations in a neighborhood of a current estimate of a rotation should be reachable by parametrization R [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. Since SO(3) is a 3D manifold, such a neighborhood is 3D and jacobian ∂ R of R should have rank 3. Thus a minimal parametrization R meets k = 3, i.e. Ω ⊆ R 3 . Unfortunately, all 3D parametrizations of SO(3) have singularities [START_REF] Singla | How to avoid singularity when using Euler angles[END_REF].

We detail the case of Euler parametrization E(α, β , γ) = R z (γ)R y (β )R x (α) where R x (α), R y (β ) and R z (γ) are the rotations about axes x-y-z with angles α-β -γ. The singularities of E are the points (α, β , γ) ∈ R 3 where ∂ E is rank deficient, i.e. planes β = π/2 + pπ such that p ∈ Z (a detailed proof is in the supplementary material). Parametrization E can be used in BA if β is as far as possible to π/2 + πZ, e.g. local Euler angles [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] such that |β | 1.

The angle-axis parametrization used in [START_REF] Oth | Rolling shutter camera calibration[END_REF] has a bounded angle domain due to singularities at angles 2πZ * . This restricts the continuous camera motion: rotations around a fixed axis must have angle range in ] -2π, 2π[, e.g. no more than two full turns around a building.

Keep Away from the Singularities of Euler Parametrization

Here we define R using E and keep away from all singularities. According to Sec. 3.2, R is a global parametrization for all rotations of the continuous camera motion M(t). Note that there is a single R for all times, in contrast to the local methods in [START_REF] Singla | How to avoid singularity when using Euler angles[END_REF][START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] which switch R over time. Furthermore, R has singularities since it is minimal (Sec. 4.1). Thus we do an assumption on the camera motion to keep away from the singularities. Now we remind that the multi-camera is helmet-held. All yaw motions of the head are possible since the user can move in all horizontal directions. We assume that the pitch and roll of the head are small, i.e. the viewing direction of the user is roughly pointing toward the horizon without odd roll rotations. We believe that this assumption is reasonable for an user exploring the environment without a special objective like grasping at objects on the ground.

Let R 0 i be the rotation of the initial m i computed by standard BA. Thus all (R 0 i ) T R 0 j are roughly rotations sharing a same axis v ∈ R 3 . Let R(v, θ ) be the rotation with axis v and angle θ . There are rotation R and angles θ i such that R 0 i ≈ RR(v, θ i ). Let k = 0 0 1 T . Let rotations A and B be such that ∀i,

AR 0 i B ≈ R(k, γ i ). Angles (α i , β i , γ i ) meet E(α i , β i , γ i ) = AR 0
i B and 2π multiples can be added to them. We would like that α i , β i , γ i are coordinates of M(t i ). We choose β i that has the smallest

|β i |. Since E(α i , β i , γ i ) ≈ R(k, γ i ), |β i |
is small enough to keep away from E singularities. We also remind that M is continuous 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 (|β i -β i-1 | is also small, and we do similarly for α i ). Last we define R(α,

β , γ) = A -1 E(α, β , γ)B -1 to obtain R(α i , β i , γ i ) = R 0 i .
Since R has the same singularities than E (supplementary material) and β ≈ β i during our BA, (α, β , γ) is far from the R singularities.

Technical Details: Estimate A and B

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

i ≈ R(v, θ i ). Since R 0 i ≈ RR(v, θ i ), R 0 i v ≈ Rv. Let ṽ = ∑ i R 0 i v/|| ∑ i R 0 i v||. Thus ṽ ≈ Rv ≈ R 0 i v. Let R be a rotation such that Rṽ = v. Since RR 0 i v ≈ Rṽ = v, RR 0 i ≈ R(v, θ i ). Let R be a rotation such that R v = k. We obtain R RR 0 i R T ≈ R(k, γ i ). Thus A = R R and B = R T .
Figure 1: Two multi-cameras formed by four GoPro Hero3 cameras and images taken at a viewpoint for every dataset. Top and middle: the cameras are enclosed in a cardboard (for small baseline). Bottom: the housings provided with the cameras are used.

Experiments

Datasets and Notations

The multi-camera is defined by four GoPro Hero3 cameras (Fig. 1) that have the same setting at a time, except the camera gain that evolves independently for every camera. We assume that time offsets and calibrations do not change in a video. Tab. 1 summarizes our datasets: three real (multi-camera) videos under various conditions (bike riding in a city, walking in a town [START_REF] Lhuillier | Synchronization and self-calibration for helmet-held consumer cameras, applications to immersive 3d modeling and 360 videos[END_REF], paragliding flying at very low height above a hill) and one synthetic video. In all cases, a 360 • field of view around the head is obtained. The ground truth of line delay τ is available thanks to a strobe. Furthermore, we use global shutter and central approximations [START_REF] Lhuillier | Synchronization and self-calibration for helmet-held consumer cameras, applications to immersive 3d modeling and 360 videos[END_REF] to obtain the initial self-calibration (FA synchronization, intrinsic parameters, 3D points, multi-camera poses and relative poses) and a concurrent SFA synchronization based on instantaneous angular velocity (we call it Sync).

BikeCity2 is generated by ray-tracing of a synthetic urban scene having real textures and by moving the camera along a trajectory that mimics that of BikeCity1 (the "pose noise", i.e. relative poses 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". BikeCity2 has ground truth: f ∆ 0 = 0, f ∆ 1 = 0.25, f ∆ 2 = 0.5, and f ∆ 3 = 0.75 where f = f ps (if f ∆ j = 1, ∆ j is the time between two consecutive frames).

We note that the setting of the cameras of FlyHill is different to those of the other videos (frequency, resolution, orientations, baseline). The baseline in the others is as small as possible for the central approximation. Furthermore, FlyHill is the most difficult case due to rolling shutter. Indeed, its fps is twice smaller and it includes faster head turns.

We use notations: C (central approximation) estimates R j and fixes t j = 0, NC (noncentral) estimates both R j and t j , RS estimates τ, SFA estimates ∆ j , GS means τ = 0 and FA means ∆ j = 0. For example, a GS+SFA+NC bundle adjustment fixes τ = 0 and estimates 1: Our videos: angular resolution r (milliradians), diameter b of multi-camera centers, line delay τ (ground truth), trajectory length l, numbers of frames f r and keyframes k f r and 3D points #3D, maximum of angles |β i | (radians) of our parametrization in Sec. 4.2. simultaneously all ∆ j , R j , t j , the keyframe poses m i and the 3D points. The threshold for inlier selections is set to 4 pixels in all videos.

Name f ps r (mr) b (cm) τ (µs) l (m) f r k f r #3D ||β i || ∞ BikeCity1

One New Assumption at Once

There are several new assumptions NC, SFA and RS. Thus we first examine what they provide separately: we experiment GS+NC+FA, GS+C+SFA and RS+C+FA bundle adjustments. The number of independent and optimized parameters of these BAs are x + 9, x + 3 and x + 1 (respectively), where x is the number of parameters of the initial GS+C+FA BA.

In Tab. 2, the inlier set is fixed in every video to compare the improvements in term of RMS of reprojection errors (in pixels). The RMS decreases are small (less than 1.7%), except for FlyHill: GS+C+SFA has 3.9% and RS+C+FA has 3.2%. This confirms that the RS and SFA effects are non negligible for FlyHill due to the fast image motion (it is faster than in other videos). We also see that the NC assumption has the lowest impact on the RMS in spite of its larger number of parameters. At this point, the relative error of τ is quite large for FlyHill (58%), it also important for the others (3.8%-13.3%). The error of f ∆ j (BikeCity2), or the difference between our f ∆ j and those of Sync (others), can reach 0.14 in BikeCity2 or 0.24 in WalkTown. Such discrepancies look large since we expect that | f ∆ j | ∈ [0, 1], but the resulting discrepancies for the 3D locations of the multi-camera are small. For example, the mean distance between multi-camera poses for consecutive images of WalkTown is 900/70300m, thus f ∆ j = 0.24 implies a 3D discrepancy of only 2.9mm.

We continue these experiments by alternating inlier updates and BAs in Tab. 3. Then τ is improved: the relative error is less than 7.9% except for FlyHill (42%). The f ∆ j estimation is similar to that in Tab. 2 (the RMS of differences of f ∆ j is less than 0.05). The inlier sets increases slightly: 0.6-0.9% for RS+C+FA and GS+C+SFA of FlyHill and less than 0.1% elsewhere; their main computations have been done before by initialization BA (GS+C+FA).

Several New Assumptions at Once

Now we try RS+SFA simultaneously and study differences between these results and the previous ones in Tab. 3 and Sync. More precisely, we compute RS+NC+SFA and RS+C+SFA like this. Get the GS+NC+FA result in Tab. 3, apply RS+NC+FA and then RS+NC+SFA. Get the GS+C+SFA result in Tab. 3, apply RS+C+SFA (all BAs include inlier updates). Fig. 2 shows views of the RS+NC+SFA results.

We first focus on BikeCity2 (Tab. 4). The accuracy of ∆ j is better than that in Tab. 3 and Sync (the RMS of f ∆ j difference between Tab. 4 and ground truth is 0.027 for RS+NC+SFA or 0.019 for RS+C+SFA). However the relative error of τ is worse: 6.8%-8.5% (0.31% in Tab. 3). Thus the simultaneous use of RS+SFA improves ∆ j but provides a bias for τ.

Method i-RMS f-RMS f ∆ 1 f ∆ 2 f ∆ 3 10000 f τ BikeCity1 Sync or G.T.
- 3: Results of our BAs with increasing set of inliers for every video (#3D is the number of 3D inlier points added to those in Tab 2, i.e. Tab 1). Percentages are relative errors.

For FlyHill, the relative error of τ is less that 5.1% and is quite better than that in Tab. 3; offsets ∆ j of RS+NC+SFA and RS+C+SFA are similar. Once more, the largest increase of inlier set (3.5%) is observed in this video. The relative error of τ is better for BikeCity1 (less than 5.3%), but it is worse for WalkTown (less than 9.8%). For FlyHill and BikeCity1-2, the RMS of f ∆ j diff. between Tab. 4 and Sync (range from 0.1 to 0.2) is greater than the RMS of f ∆ j diff. between Tab. 3 and Sync (range from 0.05 to 0.08). Both Sync and GS+C+SFA (Tab. 3) use the same GS+C assumption and we believe that they provide similar ∆ j for this reason (this does not mean that the ∆ j in Tab. 3 are better than the ∆ j in Tab. 4).

Method

#3D f ∆ 1 f ∆ 2 f ∆ 3 10000 f τ BikeCity1 RS+NC+SFA +492 -0.

Parametrization of the Multi-Camera Orientations

Here we examine the Euler angles involved in our rotation parametrization in Sec. 4.2. Fig. 3 illustrates function E that maps keyframe number i to (α i , β i , γ i ) for BikeCity1. Function E looks continuous (zoom in to see the blue crosses); the largest value of |γ iγ i-1 | is equal to 0.61 rad. Such a result is expected since E M is assumed to be C 3 continuous and the keyframe sampling t i is dense enough to obtain a successful SfM result. Here t i+1 -t i ranges from 0.1s to 2.3s. Furthermore, |β i | is as small as to keep from the singularities β ∈ π/2 + πZ. According to Tab. 1, all |β i | are less than 0.44 rad. The |β i | RMS is about 0.3-0.6 times the |β i | maximum for every video.

Last we detail consequences of a naive use of Euler angles ignoring singularities (using notations in Sec. 4). Assume that the initial multi-camera poses meet R 0 i ≈ R z (γ i )R y (π/2). This is possible because of a "bad" choice of coordinate systems: if the world coordinate system is rotated by A and the multi-camera coordinate system is rotated by BR y (π/2), R 0 i is replaced by AR 0 i BR y (π/2) ≈ R z (γ i )R y (π/2). Now we naively set R = E and redo the experiments of BikeCity1. Then β i ≈ π/2 (this is very close to a singularity). Angles α i and γ i are quite more perturbed although they are chosen such that function E is as continuous as possible (max i |γ iγ i-1 | = 3.11 and max i |α iα i-1 | = 3.13). The new values of 10000 f τ

Conclusion

We present the first bundle adjustment for multi-cameras that estimate not only rolling shutter (line delay) but also synchronization (time offsets), in addition to the usual 3D parameters (points, camera and multi-camera poses). In contrast to the previous Structure-from-Motion methods involving rolling shutter, only keyframes are involved and we deal with larger trajectories (600m-2.5km). The multi-camera motion is modeled at all times thanks to Taylor approximations and a careful use of Euler angles avoiding singularities. We experiment in cases that we believe useful: several and identical consumer cameras mounted on a helmet.

At first glance, our approximations seem hazardous if the user does a motion that is not consistent with the neighboring keyframes. Anyway, the majority of keyframes provides accurate enough approximation to obtain the following results in our non trivial datasets. The relative error of the estimated line delay is less than 7.9% except in the most difficult case with faster head motions; the simultaneous estimation of line delay and time offsets can provide bias but it also provides the best result (5.1%) for the most difficult case. The best (subframe-accurate) time offsets are given by the simultaneous estimation.

Several extensions are possible: adding estimated parameters only for keyframes where our model of multi-camera motion is not accurate, adding parameters by taking care of overfitting (e.g. intrinsic parameters, one line delay and frame rate per camera), trying alternative camera models and rotation parametrizations, improving applications like 3D modeling and 360 video.

Figure 2 :

 2 Figure 2: From left to right: reconstructions of BikeCity1, WalkTown, FlyHill and BikeCity2 by RS+NC+SFA without loop closure. The FlyHill trajectory has a lot of sharp S turns.

Figure 3 :

 3 Figure 3: Euler angles for BikeCity1.

Table 2 :

 2 BA results for a fixed set of inliers for every video. i-RMS and f-RMS are RMS before and after BA.

					0.042 -0.163 0.306	9.120
		GS+NC+FA 0.9550 0.9534
		GS+C+SFA 0.9550 0.9520 -0.122 -0.147 0.314
		RS+C+FA	0.9550 0.9484			7.898
	WalkTown Sync or G.T.			0.517	0.474	0.443	9.120
		GS+NC+FA 0.9452 0.9406
		GS+C+SFA 0.9452 0.9391 0.715	0.714	0.510
		RS+C+FA	0.9452 0.9391			8.714
	FlyHill	Sync or G.T.			0.207	-5e-4 -0.358	5.424
		GS+NC+FA 1.3643 1.3537
		GS+C+SFA 1.3643 1.3107 0.058 -0.024 -0.231
		RS+C+FA	1.3643 1.3207			2.270
	BikeCity2	G.T.			0.25	0.5	0.75	9.120
		GS+NC+FA 0.8124 0.8121
		GS+C+SFA 0.8124 0.7985 0.388	0.476	0.788
		RS+C+FA	0.8124 0.8009			8.777
		Method	#3D	f ∆ 1	f ∆ 2	f ∆ 3	10000 f τ
	BikeCity1 GS+NC+FA	+87
		GS+C+SFA	+130 -0.137 -0.156 0.336
		RS+C+FA	+323				8.401 (7.9%)
	WalkTown GS+NC+FA +165
		GS+C+SFA	+235	0.750	0.746	0.535
		RS+C+FA	+242				9.020 (1.1%)
	FlyHill	GS+NC+FA +448
		GS+C+SFA +5231 0.089 -0.028 -0.312
		RS+C+FA	+3627				3.153 (42%)
	BikeCity2 GS+NC+FA	+9
		GS+C+SFA	+38	0.403	0.495	0.836
		RS+C+FA	+21				9.092 (0.31%)
	Table

Table 4 :

 4 Results of our BAs with increasing set of inliers for every video.

				357 -0.155 0.151	8.922 (2.1%)
		RS+C+SFA	+379	-0.355 -0.148 0.156	8.636 (5.3%)
	WalkTown RS+NC+SFA	+489	0.543	0.807	0.339 10.016 (9.8%)
		RS+C+SFA	+343	0.567	0.795	0.343	9.480 (3.5%)
	FlyHill	RS+NC+SFA +19659 0.286	0.200 -0.326 5.700 (5.1%)
		RS+C+SFA	+19608 0.287	0.200 -0.330 5.595 (3.1%)
	BikeCity2 RS+NC+SFA	+53	0.256	0.542	0.772	8.497 (6.8%)
		RS+C+SFA	+50	0.251	0.530	0.762	8.342 (8.5%)
		Sync		0.404	0.398	0.829


Vadim Litvinov

Maxime Lhuillier

Incremental Solid Modeling from Sparse Structure-from-Motion Data with Improved Visual Artifacts Removal

Keywords: 3D shape recovery, Stereo and multiple view geometry, Reconstruction and camera motion estimation

In the recent years, a family of 2-manifold surface reconstruction methods from a sparse Structure-from-Motion points cloud based on 3D Delaunay triangulation was developed. This family consists of batch and incremental variations which include a step that remove visual artifacts. Although been necessary in the term of surface quality, this step is slow compared to the other parts of the algorithm and is not well suited to be used in an incremental manner. In this paper, we present two other methods for removing visual artifacts. They are evaluated and compared to the previous one in the incremental context where the need of new methods is the highest. Taken separately, they provide medium results, but used together they are as good as the old method in the terms of surface quality, and at the same time, processing time is almost three times smaller.

I. INTRODUCTION

Despite being around for quite a long time, surface reconstruction is still an active research topic. In the last few years, several methods were developed to perform solid modeling from a sparse Structure-from-Motion (SfM) point cloud, both in batch [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF] and incremental [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] frameworks. The main interest of this family of methods is their time and space complexities due to the fact that a dense stereo step is unnecessary. Furthermore, this series of algorithms produces a 2-manifold output which is useful to initialize a more precise, but heavy, dense reconstructions or simply allow to perform efficient surface smoothing. Other surface reconstruction methods [START_REF] Pan | ProFROMA: Probabilistic feature-based on-line rapid model acquisition[END_REF]- [START_REF] Hoppe | Incremental surface extraction from sparse Structure-from-Motion point clouds[END_REF] work in a similar (sparse) context but produce non manifold surfaces.

For memory, a surface has the 2-manifold property if and only if every point of it has a neighborhood which is topologically a disk, i.e. each triangle of the surface is exactly connected by its three edges to three other triangles, the surface has neither holes nor self-intersections and it divides the space into separate regions. If the surface is 2-manifold, its genus is defined as the number of its handles. For example, the genus of a sphere is 0, that of a torus is 1, etc. The incremental method [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] considered in this paper is a sculpting method in a 3D Delaunay triangulation. It takes a sparse cloud of 3D points sampled on an unknown surface along with their visibility information produced by SfM. The 3D Delaunay triangulation T of P is a partition of the convex hull of P such as the circumscribing sphere of every tetrahedron does not contain a vertex in its interior. All the tetrahedra of T are labeled free-space or matter using the visibility Fig. 1. A visual artifact (left) and its removal (middle and right). The triangle orientations are colored (white: ground, red-green-blue: vertical, black: sky). constraints. A ray is defined as a line segment between a camera and a point observed by this camera. A tetrahedron is free-space if it is intersected by at last one ray.

The list of triangles separating free-space and matter can, at first glance, be directly considered as the output surface, but the resulting surface is not guaranteed to be a manifold. To enforce the topology constraint, a region growing approach as in [START_REF] Faugeras | Representing stereo data with the Delaunay triangulation[END_REF] can be used, but alone it is limited to genus 0 surfaces. So methods [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF], [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] expand it by a topology extension step.

However, the region growing and topology extension steps alone suffer from a visual artifacts problem as illustrated on Fig. 1. The artifact that connects the wall of the building to the ground on the left part of the figure does not exist in reality and should be removed as shown on the central part. This is only one example of situations when the problem arise. It can also occur in other contexts [START_REF] Chaine | A geometric convection approach of 3D reconstruction[END_REF].

There were several previous attempts to solve this problem. First of all, computational geometry only approaches [START_REF] Wood | Removing excess topology from isosurfaces[END_REF], [START_REF] Zhou | Topology repair of solid models using skeletons[END_REF] that use scanner data consider that spurious handles (a particular, but very usual kind of artifact) are usually small by contrast to the real world handles that would be large. Unfortunately this is not always true in our case. Another idea [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF] is to use the visibility information provided by SfM.

The artifact removal method in Sec. 4 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF] is efficient in terms of resulting surface quality, but, unfortunately, it is the slowest part of the calculation if the input is large enough. Another limitation is the fact that this algorithm has an artificial (Steiner) points addition step that breaks Delaunay property of the triangulation. Although this is not really a problem in the batch case [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF], this is not recommended in the incremental framework [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] where the triangulation is updated by vertex additions after every artifact removal step.

Indeed, the Delaunay triangulation has theoretical warranties for surface reconstruction in both geometric [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] and topological [START_REF] Amenta | Surface reconstruction by Voronoi filtering[END_REF] senses. We also benefit by the fast and incremental 3D Delaunay triangulation in the CGAL library [12].

The main contribution of this paper is two other visual artifacts removal methods that are faster than the previous one and don't require Steiner points addition. Another contribution is a comparison between these methods applied in the incremental [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] framework.

Sec. II summarizes the batch and incremental surface reconstruction methods. We remind previous visual artifacts removal method in Sec. III-A and then explain in detail the two new methods proposed by this paper in the remaining of Sec. III. Finally, we compare these methods in Sec. IV.

II. SUMMARY OF SOLID MODELING METHODS

In this section, we summarize the batch and incremental solid modeling methods. The first and common step of the two algorithms is Harris points extraction and their 3D reconstruction by a Structure-from-Motion algorithm. Let P be a set of 3D points reconstructed by SfM, C a list of camera locations. Let R be the list of rays defined as follow: ∀c i ∈ C, ∀p j ∈ P , segment c i p j ∈ R if and only if c i has observed p j .

A. Batch surface reconstruction method

Here is a brief overview of the batch method in [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF]. First of all, a 3D Delaunay triangulation T of P is constructed (T is a list of tetrahedra). For all ∆ ∈ T , we define I(∆) as the number of intersections between ∆ and the segments of R. Each tetrahedron ∆ ∈ T is labeled free-space or matter: ∆ is free-space if and only if I(∆) = 0. We compute the list F ⊆ T of the free-space tetrahedra.

As the next step, we define another partition of T . Let O ⊆ F such as the border ∂O of O is 2-manifold (∂O is the list of triangles which are faces of exactly one tetrahedron in O). The tetrahedra in O are outside and the others are inside.

Initially O = ∅. We begin by the region growing step: for each ∆ ∈ F , we add ∆ to O if ∂O remains manifold. This property is checked using a very fast test. We begin by adding tetrahedra with highest I(∆). This way, we begin by adding tetrahedra with a highest chance of being free-space.

Then, we use the topology extension step to allow the genus of the resulting surface be higher than zero. We add the tetrahedra ∆ ∈ F by packs of tetrahedra incident to a common vertex v ∈ P . To ensure the manifold property, in this configuration, a slower but more general test is used.

We repeat these two steps until no more tetrahedra can be added to O. The resulting surface is the border ∂O of O. Some post-processing steps are then applied to achieve better surface quality (more details in [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF]), including the artifact removal method which is summarized in Sec. III-A for paper clarity.

B. Incremental surface reconstruction method

To experiment our artifacts removal methods in an incremental context, we use the algorithm described in [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF]. Let P t be the set of points computed at time t that will not be modified by SfM after t (p ∈ P t and c i p ∈ R imply i ≤ t; i = t exists). At the keyframe/time t + 1, the algorithm has the following input: T t is a 3D Delaunay triangulation of P t , F t ⊆ T t is a list of free-space tetrahedra of T t and O t ⊆ F t is a list of outside tetrahedra such that ∂O t is a 2-manifold.

The goal of the algorithm is to update the surface with a set of new SfM points P t+1 . They can not directly be added into the triangulation T t because this can modify O t and break the manifold property of its border. So we compute a ball B centered on c t+1 with the radius big enough to enclose P t+1 and neighboring tetrahedra (see Fig. 2). Then we remove tetrahedra from O t one-by-one and by packs in the similar manner than Sec. II-A until ∂O t is still a manifold and we get as close as possible to O t ∩ B = ∅. This is the shrinking step.

As the next step, we safely add P t+1 into T t without modifying O t . We obtain T t+1 . Then, we update the freespace/matter status of tetrahedra and obtain F t+1 ⊆ T t+1 .

Finally, to compute the new surface, we initialize O t+1 = O t and we apply the region growing and topology extension algorithms described in Sec. II-A. To enhance the resulting surface quality, some post-processing steps are applied locally to the updated region (see [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] for details).

III. VISUAL ARTIFACTS REMOVAL METHODS

Before we begin the discussion about different artifacts removal methods, we should precisely define what we seek to remove. Using the definitions in Sec. II, the 3D Delaunay triangulation T is encoded by an adjacency graph Γ T in which the nodes are the tetrahedra of T and the edges are the triangles between two tetrahedra. In the same manner, for any set S ⊆ T of tetrahedra, Γ S ⊆ Γ T is the corresponding adjacency graph.

We define a visual artifact

A by A ⊆ F \ O and Γ A is connected, i.e. a set of tetrahedra that is included in the inside volume (A ∩ O = ∅), but not in the scene matter (A ⊆ F ).
Unfortunately, detecting all the artifacts in the resulting model using this definition alone would be too slow to be useful in practice. So we define a visually critical edge [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF] as a segment ab such as a ∈ P, b ∈ P and ∃c ∈ C such as acb > α where α is a user defined threshold (α = 5 deg in this paper). Then, we define a visually critical artifact as a visual artifact which has (at least) a tetrahedron containing a visually critical edge.

The artifacts removal methods discussed in this section seek to remove as many visually critical artifacts as they can. We do not use the term spurious handle as in [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF], because it is too restrictive. Our algorithms deal with more than just "handles".

0 2 2 0 O 1 3 0 2 2 0 O 2 3 0 2 2 0 O 3 f (O 3 ) > f (O 1 ) > f (O 2 )

A. Summary of previous method

In this section, we briefly discuss the visual artifacts removal method described in [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF]. As its entry, this algorithm takes the list E α of the visually critical edges.

First of all, the method splits each edge of E α in its middle by adding a Steiner point (i.e. artificial point which wasn't seen by any camera). Each edge is split into two parts and the resulting edges are stored in a list E split . Each tetrahedron including the edge is also split in two and the two resulting tetrahedra are assigned the same number of intersections and status as the initial one.

This way, we reduce the size of tetrahedra and hope to locally unlock region growing in the neighborhood of this edge. Furthermore, ∂O is still manifold. The drawback is that the triangulation is not guaranteed to be Delaunay anymore.

After this first step, we create a list V e of the end vertices of the edges of E split . For each vertex v ∈ V e , we apply a force/repair cycle as following.

We begin by the force step. Let L v be the list of the tetrahedra incident to v. We define G = (L v ∩ F ) \ O. So G is a set of tetrahedra forming a visual artifact. We force them into outside: O ← O ∪ G. If the border of O remains manifold, we get what we wanted. If it doesn't, it contains some vertices whose neighborhood is not topologically a disk. These vertices are called singular vertices in the remaining of this paper. We call n the number of these singular vertices. So, we try to repair O using the next step.

Then, we execute the repair step. We seek to add a bunch of free-space tetrahedra to O such that ∂O become manifold once again. To achieve this goal, we apply a local region growing algorithm to O in F starting in the neighborhood of G and which decreases n. The algorithm stops if a number of iterations g 0 (fixed by the user) has been reached or no more tetrahedra can be added to O. If the repair step succeeds (i.e. n = 0: ∂O is manifold), we was able to remove a visual artifact and proceed to the next vertex. Otherwise, we restore O to the previous state.

B. New visual artifacts removal method

As was previously noted, the method in Sec. III-A is good to improve the surface quality, but unfortunately it has two problems. The first is the need of Steiner points for good results, which break the Delaunay property of the triangulation.

The second problem is that artifacts removal method is slow (48%-72% of the total processing time according to Tab. 2 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF]). So, we propose a faster visual artifacts removal method that doesn't require these Steiner points.

The region growing described in Sec. II is a greedy optimization algorithm that maximizes

f (O) = ∆∈O I(∆) (1)
in the discrete search space of tetrahedra lists O such that O ⊆ F and ∂O is a 2-manifold. In [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF], we have a steepest descent heuristic for function -f (O) (by region growing of O) which can get stuck to a local maximizer of f (O). A visual artifact (e.g. spurious handle) can be seen in this situation. The basic idea of our new method is to remove some tetrahedra from O (and so to decrease f (O)) to kick the algorithm out of its local extrema.

As in the previous section, we call E α the list of visually critical edges. Let G α ⊆ F such that every ∆ ∈ G α has an edge in E α . Then, for every vertex v of both ∂O and G α , we force neighboring tetrahedra out of O and we try local region growing beginning from neighboring tetrahedra included in G α . If the final value of f (O) is greater than the initial one, we are able to escape from a local maximum, otherwise we revert everything to the initial state and try another vertex v. See Fig. 3 for an example and the algorithm below. Once we tried all ∂O vertices, we complete the result using region growing and topology extension (Sec. II-A) restricted to the tetrahedra included in B (Sec. II-B).

The algorithm in pseudo-code is as follow:

1: function REGION GROWING(∆ 0 ) Local region growing in G α 2:
Let Q be a priority queue of tetrahedra based on I(∆)

3: PUSH(Q, ∆ 0 ) 4:
L add ← ∅ 5:

while Q = ∅ do 6: ∆ ← POP(Q) 7: if ∆ / ∈ O and ∂(O ∪ {∆}) is manifold then 8: O ← O ∪ {∆} 9:
L add ← L add ∪ {∆} for all vertex v of both ∂O and G α do 24:

N v ← all the tetrahedra incident to v 25:

L sub ← O ∩ N v 26: L seed ← (G α ∩ N v ) \ O 27:
L add ← ∅ end for 45:

until s old = s 46: end procedure

C. New spurious handle removal method

We also propose a visual artifacts removal algorithm that is a specialization of Sec. III-A method. The latter is slow mainly because it consists in many attempts to remove a small pack of tetrahedra and the majority of these attempts are unsuccessful. So, the idea is to remove bigger packs of tetrahedra by placing ourselves in a less general context. Instead of trying to remove all the visually critical artifacts, we seek to remove a particular kind of artifact: the handle. The fourth column of Fig. 6 shows an example. This is the most visible kind of artifact for the human eye.

The basic idea of the algorithm is to remove the handle (by contrast to Sec. III-A where we remove all the tetrahedra including a vertex) and then try to restore the manifold property of O by a local region growing. To do that, we begin as usual, by computing the list E α of visually critical edges.

For each edge ab ∈ E α , we check if it is contained in a handle.

We define a plane π perpendicular to ab and intersecting segment ab in some point. In practice we try several planes intersecting ab in 2a+b 3 , a+b 2 and a+2b 3 . Let L π be the list of the tetrahedra intersected by π. Let N ab be the list of the tetrahedra including edge ab. A handle H is a set of tetrahedra forming a visual artifact, so H ⊆ F \ O. With this definition, we begin to form our handle by

H ← (N ab ∩ L π ∩ F ) \ O.
Let N H be the list of the tetrahedra directly adjacent to the set H (i.e. ∀∆ ∈ N H , ∆ has a 4-neighbor in H). We iteratively grow H by performing

H ← H ∪ ((N H ∩ L π ∩ F ) \ O)
until no more tetrahedra can be added. Then we check that the final H is surrounded by O in plane π (see Fig. 4), i.e.

∀∆ ∈ (N H ∩ L π ) \ H, ∆ ∈ O.
Once we have detected a handle H, we try to remove it as in Sec. III-A. First we force H in O (i.e O ← O ∪H), then we try to restore the 2-manifold property of ∂O using the repair step initialized by the tetrahedra in (N H ∩F )\O. If we succeed, we remove an artifact, otherwise we undo everything and try another edge in E α .

IV. EXPERIMENTS

A. Dataset

To perform a quantitative evaluation of the visual artifacts removal methods discussed in this paper, we use the SfM dataset in [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF]. The video is a 2.5 km. long trajectory taken in a general urban environment by a PointGrey Ladybug omnidirectional camera. Our SfM is based on local bundle adjustment [START_REF] Mouragnon | Generic and real-time Structure-from-Motion[END_REF] and is refined by a global bundle adjustment to close the large loop (2.3 km). The SfM selects 1306 keyframes out of 7735 frames and detects 2.64M Harris interest points. 483k 3D points are reconstructed and shown in Fig. 5. Thus, we only have 193 points per meter to reconstruct the surface.

B. How to compare the artifact removal methods ?

Every removal method is integrated in the surface postprocessing step of the incremental surface reconstruction method [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] (before smoothing). To compare the methods in terms of output surface quality, we manually count the visually critical artifacts (Sec. III) remaining on the final surface. Because we seek to remove in priority the artifacts that are visually critical (and so are easily noticed by a human eye), we consider this number as a good quantitative metric. Fig. 6 shows examples of artifacts that are manually counted. Moreover, we also estimate the final number of freespace inside tetrahedra (i.e. the union of the visual artifacts) and the final value of the objective function f (Eq. 1). The latter quantifies the ability of every method to unlock region growing and topology extension steps (Sec. II-A), or in other words, the ability to escape from local extremum of f . Fig. 6. Results of artifacts removal methods at four locations (one location per column). Lines from top to bottom: textured scene, no artifacts removal, method III.A, method III.B, method III.C, method III.B&C.

C. Comparisons

We evaluate five visual artifacts removal methods: None (no removal method), III.A (the method in [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] summarized in Sec. III-A), III.B (escape from local extremum using the method in Sec. III-B), III.C (handle removal using the method in Sec. III-C), III.B&C (use III.C after III.B). The results are summarized in Tab. I. The removal methods III.A, III.B, III.C and III.B&C provide similar improvements (increases) of the objective function f . The differences between them are small (less than 0.034%), but we see that III.A is slightly better than III.B and III.C taken separately. Furthermore, we see that the combination of III.B and III.C is slightly better than III.A, and this is confirmed by both the number of tetrahedra in F \ O (union of all visual artifacts) and the number of artifacts that are manually detected. Fig. 6 compares the results of all methods at four locations in the reconstruction and is consistent with these comparisons. We note that III.B has important visual artifacts (as the one in the left of this figure) although it has a good (small) F \ O, and III.C can correct them in spite of its greater F \ O. Then we think that the combination III.B&C is a good choice.

A difference between the four methods is the calculation time. Indeed, III.B, III.C and III.B&C are significantly faster than III.A since their mean time per keyframe is about 2.5-4 times smaller (Tab. I and Fig. 7). We also see that III.A needs the largest number of tests that check the manifold property at a vertex.

D. More details on III.B&C Result

Fig. 8 shows a global view of the 3D models reconstructed by the method in [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF] including the artifact removal method III.B&C. At the end of the process, the surface has 600k triangles and the 3D Delaunay has 2.09M tetrahedra. The computation times for every keyframe is in Fig. 9. We use a 4xIntel Xeon W3530 at 2.8Ghz. The surface reconstruction time for every keyframe has mean 1.87 s. and maximal value 5.39 s. The largest times are at the sequence end, where the large loop is closed. As in [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF], this is due to the fact that we reconstruct new points of the loop end at a location where points are already reconstructed at the loop beginning. The joint video shows the progressive surface reconstruction and the final surface. We observe differences with [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF]: the artifact removal step (previously named handle removal) is not the most costly step anymore, and the global times are roughly divided by 1.5-2. Fig. 9. Computation times at every keyframe for the steps of our incremental surface reconstruction method including III.B&C as visual artifact removal. Note that these times are accumulated for steps shrink, . . . , artifact removal.

These improvements are the consequences of both our artifact removals (our contribution) and reasons in the Appendix.

V. CONCLUSION

This paper improves the visual artifact removal step of a recent incremental method which reconstructs a 2-manifold surface from sparse Structure-from-Motion data. We replace this step by the same method restricted on particular cases of visual artifacts (handles) and preceded by another method which escapes from the local extrema of the objective function. Then the processing time of the step is greatly reduced without loss of surface quality. We experiment on a 2.5km. long video sequence taken by an omnidirectional camera moving in an urban environment.

Future work includes improvements of calculation time and surface quality, e.g. reconstruct and integrate the image contours in the 3D Delaunay triangulation, improve the interest point matching, investigate other objective cost functions.

Fig. 2 .

 2 Fig.2. Incremental reconstruction method overview. White tetrahedra are free-space, gray tetrahedra are matter, green dots are newly added points and green dashed circle is B.

Fig. 3 .

 3 Fig. 3. Escape from local extremum thanks to our new artifacts removal (2D case). White tetrahedra are outside, gray tetrahedra are inside. The number in tetrahedron ∆ is I(∆). The green dot is the vertex considered by the algorithm and thick red lines are critical edges. Left: O before removal. Middle: force neighboring tetrahedra out of O. Right: local growing of O.

1 2 3 Fig. 4 .R

 134 Fig. 4. An example of application of spurious handle removal (in the 2D case). White tetrahedra are outside, gray tetrahedra are inside, light gray tetrahedra are freespace inside. The thick blue line is a visually critical edge, the dashed red line is plane π and red dashed triangles are selected tetrahedra.

Fig. 5 .

 5 Fig. 5. Global overview of the SfM points cloud of our test trajectory as well as an example of the Ladybug camera image.

Fig. 8 .

 8 Fig. 8. Global overview of the final surface obtained by our incremental surface reconstruction method including III.B&C as visual artifact removal. We show a textured view (top) and color coded normals (bottom). The triangles on the sky are removed for visualization purposes.

TABLE I

 I Fig. 7. Computation times at every keyframe for both previous (III.A, red) and new (III.B&C, green) visual artifact removal methods.

					.	ARTIFACTS REMOVAL RESULTS
		Removal	Mean	Max.		Num.	Size of	f=	I(∆)	Num.
									∆∈O
		method	time	time		artif.	F \ O	(M = 10 6 )	tests
			None	0	0		20	139448	29.822M	0
			III.A	1.1 s.	3.87 s	16	93259	30.012M	468M
			III.B	0.33 s.	0.94 s	18	92604	30.009M	186M
			III.C	0.27 s.	0.99 s	15	94469	30.002M	136M
		III.B&C	0.43 s.	1.2 s		13	92280	30.013M	323M
							Computation times comparison
		4								III.A III.B & III.C
		3.5
		3
		2.5
	Time (s)	2
		1.5
		1
		0.5
		0	0	200	400		600	800	1000	1200	1400
							Key frame

ACKNOWLEDGEMENTS

Thanks to CNRS, FEDER and Auvergne Région for funding.

APPENDIX

Corrections are applied to the implementation in [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF]. An implementation bug is corrected and the shrinking step is now slower but more accurate than before. We remember that the ideal result of shrinking meet O ∩ B = ∅ where B is a ball where almost all computations are done. Now this condition is not meet in only one of the 1306 keyframes (before: 4.5% of the keyframes). As a consequence, the step of the regular Steiner points grid can be reduced to 10 times mean distance between consecutive camera locations (before: 15 times). We remember that this grid is useful to reduce the size of B and these Steiner points do not break Delaunay property [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional Structure-from-Motion data[END_REF].

Moreover, the topology extension step of shrinking and growing is re-implemented in a multi-threaded fashion. As was briefly explained in Sec. II, this step consists mainly in trying to add or remove a pack of tetrahedra to the outside region. In the new implementation, we simply perform many tests simultaneously, so these steps are greatly accelerated.


Vadim Litvinov
email: vadim.litvinov@univ-bpclermont.fr

Maxime Lhuillier

Incremental Solid Modeling from Sparse and Omnidirectional Structure-from-Motion Data

This paper introduces a sparse and incremental 2-manifold surface reconstruction method. It uses a sparse 3D point cloud generated by a Structure-from-Motion algorithm (SfM) as its main input as opposed to the more common dense algorithms. Furthermore, our method is incremental: the surface is updated for every new camera pose computed by SfM, and the update occurs in a small neighborhood of the new camera pose. Compared to the other surface reconstruction methods, our method has the advantage to have all these properties at the same time. The quality and execution time of the proposed algorithm is evaluated on a large scale (2.5 km.) real sequence taken in an urban environment, and the method is quantitatively evaluated on a synthetic urban scene.

Introduction

The majority of methods which reconstruct the surface of an environment from images have a dense information as input (dense 3D point cloud or almost every image pixels). In contrast to these dense methods, the direct estimation of a surface from the sparse point cloud estimated by Structure-from-Motion step (SfM) is under-explored. The main argument against such a sparse method is that the lack of points implies an inaccurate scene surface. However, a sparse method has advantages. First, the quality (accuracy, robustness) of a 3D point is expected to be better than that of a dense stereo method, thanks to the SfM machinery involving interest point detection and bundle adjustment. Second, the resulting simplified cloud provides a simplified surface at low complexity that applications could accept: (pre)visualization with a small hardware, robot localization, initialization of a dense stereo method.

In this paper, we introduce a sparse and incremental method which reconstructs a triangulated manifold surface from SfM data. The method is incremental since the surface is locally updated for every new camera pose (and its 3D points) estimated by SfM. This is interesting for applications which require a surface while reading the video sequence. A triangulated manifold surface (shortened as "manifold") is a list of triangles in 3D such that the neighborhood of every surface point is topologically a disk. This property is needed to define surface normal and curvature [START_REF] Botsch | Polygon mesh processing[END_REF], then it is used by a lot of algorithms like surface refinement involving regularization (smoothing [START_REF] Gueziec | Cutting and stitching: converting sets of polygons to manifold surfaces[END_REF], dense stereo [START_REF] Hernandez Esteban | Silhouette and stereo fusion 3D object modeling[END_REF], . . . ) and others [START_REF] Botsch | Polygon mesh processing[END_REF][START_REF] Hoffmann | Geometric and solid modeling: an introduction[END_REF].

The majority of sparse methods are based on sculpting in a 3D Delaunay triangulation. Sec. 2 presents and discusses these methods. Then our method is described in Sec. 3 (prerequisites) and 4 (algorithm). Last Sec. [START_REF] Faugeras | Representing stereo data with the Delaunay triangulation[END_REF] provides experiments and we conclude in Sec. 6.

Sculpting Methods in a 3D Delaunay Triangulation

Let P be a set of 3D points sampled on an unknown surface. The 3D Delaunay triangulation T of P is a list of tetrahedra such that (1) the tetrahedra partition the convex hull of P, (2) their vertex set is P, (3) the circumscribing sphere of every tetrahedron does not contain a vertex in its interior. In our paper, a vertex/edge/triangle is a face of a T tetrahedron, and we use notation |L| for the union of triangles (or tetrahedra) in list L.

In the sculpting methods, a list V of tetrahedra (in T ) represents the reconstructed object whose volume is |V |. Border δV is the list of triangles which are included in exactly one tetrahedra of V . Then |δV | should be a manifold and the target surface. In the early work [START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF], P has no bad point and we should have P ⊂ |V |. This method initializes V = T , then it selects and removes a tetrahedron ∆ from V while the following conditions are met: ∆ has a triangle in δV , |δ (V \ {∆})| is manifold, P ⊂ |V \ {∆}|. An efficient rule checks that |δ (V \ {∆})| is manifold. However, the genus of |δV | is always zero, i.e. |δV | is always homeomorphic to a 2-sphere without handle/hole. The genus can be greater than 0 in [START_REF] Floriani | Managing the level of detail in 3D shape reconstruction and presentation[END_REF][START_REF] Gezahegne | Surface reconstruction with constrained sculpting[END_REF] since several tetrahedra are removed at once from V .

In our Computer Vision context, SfM provides additional visibility knowledge R i : every point p i ∈ P is computed from camera locations c j where j ∈ R i . This implies that |δV | should not intersect the rays (line segments) c j p i , j ∈ R i except at p i . The tetrahedra intersected by a ray are labeled "freespace", the others are "matter". Then the triangles between freespace and matter are good candidates to be in δV . Thanks to the visibility, we expect to reconstruct surfaces with fewer points than the methods referenced above. However, some bad points can occur since they are estimated from images, and we should not enforce P ⊂ |V |.

Several sculpting methods [START_REF] Faugeras | Representing stereo data with the Delaunay triangulation[END_REF][START_REF] Labatut | Efficient multiview reconstruction of large-scale scenes using interest points, Delaunay triangulation, and graph cut[END_REF][START_REF] Lovi | Incremental free-space carving for real-time 3D reconstruction[END_REF][START_REF] Pan | ProFORMA: probabilistic feature-based on-line rapid model acquisition[END_REF][START_REF] Yu | Incremental reconstruction of a manifold surface from sparse visual mapping[END_REF] use visibility knowledge and sparse point cloud P estimated from images. The manifold constraint on |δV | is not enforced by [START_REF] Labatut | Efficient multiview reconstruction of large-scale scenes using interest points, Delaunay triangulation, and graph cut[END_REF][START_REF] Lovi | Incremental free-space carving for real-time 3D reconstruction[END_REF][START_REF] Pan | ProFORMA: probabilistic feature-based on-line rapid model acquisition[END_REF]. Method [START_REF] Lovi | Incremental free-space carving for real-time 3D reconstruction[END_REF] is incremental, but it defines directly V as the list of the matter tetrahedra, then the resulting surface |δV | is not manifold. Method [START_REF] Pan | ProFORMA: probabilistic feature-based on-line rapid model acquisition[END_REF] is real-time (for small objects) and deals with image noise, but it is not incremental. In [START_REF] Labatut | Efficient multiview reconstruction of large-scale scenes using interest points, Delaunay triangulation, and graph cut[END_REF], a non-incremental method estimates a surface minimizing a cost involving visibility and photoconsistency. The manifold constraint is enforced in [START_REF] Faugeras | Representing stereo data with the Delaunay triangulation[END_REF][START_REF] Yu | Incremental reconstruction of a manifold surface from sparse visual mapping[END_REF]. Method [START_REF] Faugeras | Representing stereo data with the Delaunay triangulation[END_REF] is limited to genus zero surface. Method [START_REF] Yu | Incremental reconstruction of a manifold surface from sparse visual mapping[END_REF] is incremental but the time complexity of one iteration can be too large: if the camera trajectory is a loop, this complexity is at least linear to the number of camera poses in the loop.

Our contribution is an incremental sculpting method which estimates a manifold from sparse SfM data, without genus restriction and without prohibitive complexity in presence of loop in the camera trajectory. Note that other sparse methods [START_REF] Hilton | Scene modelling from sparse 3D data[END_REF][START_REF] Morris | Image-consistent surface triangulation[END_REF][START_REF] Taylor | Surface reconstruction from feature based stereo[END_REF] exist, but they can not be classified in this Section and do not estimate a manifold in an incremental scheme.

Prerequisites

Manifold Tests

According to p. 723 of [START_REF]Handbook of discrete and computational geometry[END_REF], |δV | is manifold if for every vertex v of every δV triangle, the general test is successful: the triangles in δV including v can be ordered as that t i ∩ t (i+1) mod k is an edge, and such an edge is included in exactly two triangles t i and t j . In other words [START_REF] Giblin | Graphs, surfaces and homology[END_REF], the graph of the v-opposite edges in the δV triangles must be a cycle.

t 0 , • • •t k-1 such
We complete the rule of [START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF] in our case where constraint P ⊂ |V | is not used and obtain the subtraction test. Assume that ∆ ∈ V , V ⊆ T , and |δV | is manifold. Then |δ (V \ {∆})| is manifold if the numbers of (vertices,edges,triangles) in ∆ ∩ |δV | is in {(0, 0, 0), (3, 3, 1), (4, 5, 2), (4, 6, 3), (4, 6, 4)} (Fig. 1). Only cases (3,3,1) and (4,5,2) are in [START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF].

We also need the addition test, which is similar to the subtraction test (adding ∆ to V is like subtracting ∆ from T \V ). Assume that ∆ ∈ T \V , V ⊆ T , and |δV | is manifold. Then surface |δ (V ∪ {∆})| is manifold if the numbers of (vertices,edges,triangles) in ∆ ∩ |δV | is in {(0, 0, 0), (3, 3, 1), (4, 5, 2), (4, 6, 3), (4, 6, 4)} (Fig. 1).

If we would like to add a single tetrahedron ∆ to V such that |δV | is manifold, we can use the addition test or use the general test on triangles δ (V ∪ {∆}) at every vertex in ∆ ∩ |δV |. In this case, the addition test is faster than using the general test. These tests can be converted to tests that are more convenient if T is implemented as the adjacency graph of the tetrahedra (e.g. see Appendix of [START_REF] Yu | Incremental reconstruction of a manifold surface from sparse visual mapping[END_REF] for the addition test).

Notations

At image (or time) t + 1, our method has the following input:

• 3D Delaunay triangulation T t ; • list F t of freespace tetrahedra such that F t ⊆ T t ; • list O t such that O t ⊆ F t and |δ O t | is manifold;
• list P t+1 of new Structure-from-Motion points, camera locations c t where t ≤ t + 1; and output:

• 3D Delaunay triangulation T t+1 by adding P t+1 in T t ;

• list F t+1 of freespace tetrahedra (by raytracing in T t+1 ) such that F t+1 ⊆ T t+1 ;

• list O t+1 such that |δ O t+1 | is manifold and O t+1 is the largest as possible in F t+1 .
Notations T, O, F are these tetrahedra lists if t does not need to be mentioned. The 3D ball centered at x ∈ R 3 with radius r > 0 is B r (x). Note that the vertex set of T t is

P 0 ∪ P 1 ∪ • • • ∪ P t .

Assumptions and Consequences

The following assumptions are done for complexity reasons.

First we assume that p i ∈ P t and t ∈ R i imply t ∈ R i and t ≤ t. This means that p i is reconstructed from the t-th image and previous ones, but it is not reconstructed from next ones. This implies that ray p i c t exists if p i ∈ P t and SfM does not update a reconstructed point once it is a vertex of T .

Second we assume that the lengths of all rays p i c t are bounded by r > 0. This is almost a consequence of standard SfM point filtering, since the point uncertainty roughly increases as the square of the ray length [START_REF] Lhuillier | A generic error model and its application to automatic 3d modeling of scenes using a catadioptric camera[END_REF]. Then P t is included in B r (c t ).

Third we assume that the diameters of the tetrahedra are bounded by l > 0. This condition is met if T is initialized by a large cartesian grid of Steiner vertices (i.e. extra vertices): the diameter of the circumscribing sphere of a tetrahedron is less than the diagonal length l of the (implicit) grid voxels. Since every tetrahedron created by the addition of a point has this point as vertex [START_REF] Faugeras | Representing stereo data with the Delaunay triangulation[END_REF], the addition of P t+1 is a local update of T included in B r+l (c t+1 ).

Our Incremental Sculpting Method

Assume that we add P t+1 in the Delaunay T t as soon as the O t computation is done. This destroys a list D of tetrahedra, which can contain tetrahedra in O t (Fig. 2.f). The problem is the following: if we initialize O t+1 = O t \ D, |δ O t+1 | can be non manifold, and there is no obvious method to update O t+1 such that |δ O t+1 | becomes manifold (inspired by [START_REF] Gueziec | Cutting and stitching: converting sets of polygons to manifold surfaces[END_REF], we could add new points in the Delaunay, but new tetrahedra will be destroyed and so on).

Here is our idea. List D is computed without adding P t+1 in T t . Then we initialize

O t+1 = O t , shrink O t+1 such that |δ O t+1 |

Step 1: Enclosing

We should calculate D without updating T to ensure that |δ O t+1 | is always a manifold. In practice, the results are better if we replace D by E such that D ⊆ E ⊆ T t , then shrinking will be stopped when O t+1 ∩ E = / 0 (which implies what we would like:

O t+1 ∩ D = / 0).
We chose E such that it is the list of the tetrahedra in T t included in a ball B including |D| (Fig. 2 2.a) and we remove from O t+1 tetrahedra one-at-once (using the subtraction test) or several-at-once (using the general test). During this shrinking process, |δ O t+1 | is manifold and O t+1 ⊆ F t (Fig. 2.b). The selection of the removed tetrahedra and the stopping criteria of this process are detailed in Sec. 4.7. Actually, we can not prove that the resulting O t+1 meets exactly O t+1 ∩ D = / 0, but we expect that the number of tetrahedra in O t+1 ∩ D is very small and the next step deals with this problem.

Step 3: Adding Points without Destroying δ O

We initialize T t+1 = T t and F t+1 = F t . For every point p in P t+1 , we use the following process.

Thanks to CGAL [START_REF] Cgal | [END_REF], we calculate the list D(p) of tetrahedra of T t+1 which would be destroyed if we add p in T t+1 .

If D(p) ∩ O t+1 = / 0, we add p in T t+1 (Fig. 2.c). This does not update O t+1 , thus |δ O t+1 | is still a manifold. We also apply F t+1 ← F t+1 \ D(p) i.e. the tetrahedra created by the addition of p are labelled matter. We still have O t+1 ⊆ F t+1 . The creation date t + 1 is also given to the created tetrahedra.

If D(p) ∩ O t+1 = / 0, it is difficult to update O t+1 such that |δ O t+1
| is still manifold and T t+1 is still a 3D Delaunay triangulation. We do not add p in this case, which is rare since O t+1 ∩ D is expected to be very small (Sec. 4.2). We also apply P t+1 ← P t+1 \ {p}.

Step 4: Ray-Tracing

Ray tracing labels tetrahedra in T t+1 as freespace, i.e. it increases list F t+1 . We do not trace all rays available at time t + 1 for complexity reason, so the rays are selected before tracing.

First, we trace the rays of P t+1 , i.e. every ray p i c j , j ∈ R i of every point p i in P t+1 , since this computation was not done before. Tracing ray p i c j is a walk in the adjacency graph of the tetrahedra: the walk starts from p i (which is a vertex of T ), and we go from a tetrahedron to one of its four neighboring tetrahedra if the common (face) triangle is intersected by the ray. Every tetrahedron has an intersection counter which is incremented if it is in the walk.

Second, we collect in list R every ray of P t (where t ≤ t) which can intersect a tetrahedron of N. List N is the list of the new tetrahedra, i.e. the tetrahedra which have creation date t + 1 (Sec. 4.3). The other rays do not need to be considered. Then we trace every ray in R as explained before. Since this ray was already traced, we only increment the intersection counter for the tetrahedra in N (every ray is counted once). At the end, O t+1 ⊆ F t+1 since the ray-tracing step adds new tetrahedra to the list F t+1 (Fig. 2

.d).

Here is a method to compute R efficiently. Let X N be a bounding box of the tetrahedra in N. Let X t be a bounding box of the rays of P t . For every t such that t ≤ t and X N ∩ X t = / 0, we put in R the rays of P t which intersect X N .

Step 5: O Growing

We grow O t+1 in freespace F t+1 by adding tetrahedra one-at-once (addition test) and severalat-once (general test). As required in Sec. 3.2, O t+1 ⊆ F t+1 and |δ O t+1 | is manifold (Fig. 2.e). This step is similar to that of the batch method in [START_REF] Yu | Incremental reconstruction of a manifold surface from sparse visual mapping[END_REF]. Here we give an overview of it.

First we apply a "One-Tetrahedron-at-Once Growing". A priority queue Q stores the tetrahedra in F t+1 \ O t+1 which have a triangle in δ O t+1 (we initialize Q with a tetrahedron in F t+1 ∩ E). At each step, Q provides tetrahedron ∆ with the largest ray intersection counter. We try to add ∆ to O t+1 using the addition test. If this is successful, the tetrahedra in F t+1 \ O t+1 which are adjacent to ∆ are added to Q. The process stops when Q is empty. This growing is fast thanks to the addition test, but it can not change the |δ O| genus.

Second we apply a "Several-Tetrahedra-at-Once Growing" to allow genus changes. We find a vertex v in |δ O t+1 | ∩ |E| such that all v-incident tetrahedra are in F t+1 , and try to add to O t+1 those tetrahedra which are in F t+1 \ O t+1 using the general test. If this is successful, we try to start one-tetrahedron-at once growings from these tetrahedra. The overall process stops when we can not find a successful v.

Step 6: Post-Processing

First of all, we apply a genus refinement procedure ("Handle Removal") to O t+1 ∩ E as described in [START_REF] Yu | Genus refinement of a manifold surface reconstructed by sculpting the 3D-Delaunay triangulation of structure-from-motion points[END_REF]. The only difference is that we don't add Steiner points in the middle of the edges such that T t+1 is Delaunay. Then we apply an incremental surface smoothing step [START_REF] Yu | Incremental reconstruction of a manifold surface from sparse visual mapping[END_REF]. First we apply a "One-Tetrahedron-at-Once Shrinking". We remove from Q the tetrahedron ∆ which has the smallest intersection counter. If ∆ / ∈ O t+1 or ∆ does not have a triangle in δ O t+1 , we take another ∆ in Q. Then we try to remove ∆ from O t+1 such that |δ O t+1 | remains manifold using the subtraction test (Sec. 3.1). In case of success, we add to Q the tetrahedra of O t+1 ∩ E which are adjacent to ∆. We continue until Q = / 0. This shrinking is fast thanks to the subtraction test, but it can provide

Details on the Shrinking Step

E ∩ O t+1 = / 0 (e.g. if E ∩ O t+1 = / 0 implies that |δ O t+1 | genus changes).
Second we apply a "Several-Tetrahedra-at-Once Shrinking" to allow genus changes. We find a vertex v which is both in a triangle of δ O t+1 and in a tetrahedron of E, define L as the list of tetrahedra in O t+1 ∩ E which are incident to v, apply O t+1 ← O t+1 \ L, and apply the general test for δ O t+1 at every vertex of L. In case of success, we redefine a list Q with the adjacent tetrahedra of L, and redo the one-tetrahedron-at-once shrinking above. In case of failure, we apply O t+1 ← O t+1 ∪ L and try another v. The overall process stops when we can not find a successful v.

Experiments

Real Image Sequence

The City sequence is taken by a PointGrey Ladybug omnidirectional camera. Ladybug is a (non-central calibrated) rigid multicamera system consisting of six synchronized pinhole cameras each of which takes 1024 × 768 images at 15 fps. It is mounted on a car and is about 4 meters above the ground. The trajectory is 2.5 km long and includes a large loop (2.3 km). First we provide details on SfM. A central approximation simplifies the calculation; we enforce the centers of the 6 cameras to be the same. Moreover, the size of the input images are reduced by a half to improve the computation speed. Then we adapt the incremental SfM in [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] based on local bundle adjustment. The loop is closed thanks to a global bundle adjustment. At the end of the SfM process, 1306 keyframes are selected in the 7735 available 6-tuples of images, 483k points are reconstructed from 2.64M Harris points detected and matched in the images. Thus the surface should be reconstructed with a small number of points (about 193 per meter). Fig. 3 shows a top view of the SfM result and sequence images.

Second we explain how to initialize our incremental method. Let t be the keyframe index (same notation as in Sec. 3 and4). We construct the 3D Delaunay triangulation T 40 from the points in P t where t ≤ 40, estimate freespace F 40 by ray-tracing every ray in these P t (Sec. 4.4), and obtain O 40 by the growing step (Sec. 4.5) and handle removal [START_REF] Yu | Genus refinement of a manifold surface reconstructed by sculpting the 3D-Delaunay triangulation of structure-from-motion points[END_REF]. The points of P t are filtered as follows. Let v be the vertical direction (||v|| = 1) and s = ∑ 39 t=0 ||c t+1 -c t ||/39. Every point p i in P t meets -3s ≤ (p i -c t ).v (we remove points below the ground surface), ||p i -c t || ≤ 15s (Sec. 3.3) and p i has two rays forming an angle larger than 10 • (standard filtering using aperture angle). The step of the grid of Steiner points is 15s.

Third we apply the incremental method for 40 ≤ t ≤ 1306. Fig. 4 shows the final surface as well as local views of the surface at several times t. The final surface has 528k triangles. The joint video shows the progressive surface reconstruction and the final surface (this video is also in http://www.youtube.com/watch?v=w1AQfvhGx5I).

The execution time for every keyframe is on the left of Fig. 5. We use a 4xIntel Xeon W3530 at 2.8 GHz (multi-threading is only used by the ray-tracing step). The most costly step is the handle removal (Sec. 4.6). This step and the (main part of the) other steps are only applied in E. Remind that E is the list of the tetrahedra in the enclosing ball B, which has a bounded radius (Sec. 4.1). According to Fig. 5 (right), the size of E is less than 120k, which is quite smaller than 2066k, the number of tetrahedra in the final 3D Delaunay triangulation.

The large loop (2.3 km) is closed at the very end of the sequence, where we observe a computation time increase. This is due to point increase in B: we reconstruct new points of the loop end at a location where points are already reconstructed at the loop beginning.

According to Sec. 4.2, we can not prove that O t+1 meets O t+1 ∩ D = / 0 and we expect that O t+1 ∩ D is small enough to reject a minority of points (Sec. 4.3). In our experiments, O t+1 ∩ D = / 0 is not met in only 59 keyframes out of 1306 (4.5%). For these keyframes, the maximum percentage of rejected points is 26.6% and 56 keyframes have less than 13.3% of rejected points. Then a minor amount of points can not be added in T t+1 due to O t+1 ∩ D = / 0.

Synthetic Image Sequence with Ground Truth

This Section provides the reconstruction error of our sparse scheme (central SfM + incremental surface reconstruction) applied in a synthetic urban scene, whose textures are extracted from real images taken in a city. The synthetic sequence has 1553 images generated by ray-tracing and the non-central camera in Sec. 5.1. The camera trajectory is a 621 m long and closed loop.

The SfM reconstructs 346 keyframes and 145k points from 775k Harris points detected and matched in the images. Then we apply our incremental surface reconstruction and estimate the error between the final estimated surface (182k triangles) and the ground truth surface. Now we explain how to estimate the reconstruction error. First, a pose-based registration is used to set both surfaces in the same coordinate system. The i-th keyframe has estimated location c i and ground truth location c g i ; the latter is the mean of the 6 camera centers in the ground truth coordinate system. Then we estimate the similarity transformation S minimizing E(S) = ∑ 345 i=0 ||S(c i ) -c g i || 2 , and use S to map the estimated surface in the ground truth coordinate system. Second, a ray-tracing approach is used to compute the error. Let p be a pixel in the i-th keyframe. Let p e be the first intersection of the estimated surface and the ray defined by p and the estimated i-th pose. Let p g be the first intersection of the ground truth surface and the ray defined by p and the ground truth pose of the i-th keyframe. If both p e and p g exist, we define error e(p) = ||p e -p g ||.

The pose-registration provides E(S)/346 = 14 cm. Then we uniformly sample pixels in the keyframe sequence and examine the distribution of e(p) for 769k pixels. The x-quantile q x is the real such that x percents of the e values are less than q x . We have q 10 = 12, q 20 = 16, q 30 = 18, q 40 = 22, q 50 = 28, q 60 = 48, q 70 = 76, q 80 = 112, q 90 = 174, and 3% of the e values are larger than 600 (all numbers are given in centimeters).

Conclusion

In this paper, we have introduced a new incremental 2-manifold surface reconstruction method taking a sparse SfM data as input. Compared to the previous similar algorithms, it has the advantage to produce a manifold without genus restrictions and without prohibitive complexity in case of large loops in the input trajectory. It has been experimented on a large scale real urban scene as well as on a set of synthetic data with ground truth.

Several improvements of our method are subject of future work. First, the choice of a better (smaller) enclosing area where almost all calculations are done would decrease the computation time. Second, our implementation and matching method should be improved. Moreover, several lines of investigation exists to enhance the quality of the reconstructed surface. In particular, the usage of the other type of the input primitives than the interest points should be explored. We could also think about surface denoising methods that make a better usage of the Structure-from-Motion properties and scene priors.

Figure 1 :

 1 Figure 1: Manifold tests (best viewed with colors). See the text for triplets explanation.

Figure 2 :

 2 Figure2: Sculpting method general overview (2D case). White triangles are freespace (in F), grey triangles are matter, red line is border δ O, green dots are newly inserted points and the green circle encloses the tetrahedra that will be deleted by the addition of the green dots.

  We initialize O t+1 = O t and progressively remove tetrahedra in O t+1 such that |δ O t+1 | remains manifold until O t+1 ∩ E = / 0. The O t+1 shrinking is an inverse of growing in Sec. 4.5. Let Q be the list (priority queue) of tetrahedra in O t+1 ∩ E which have a triangle in δ O t+1 .

Figure 3 :

 3 Figure 3: Camera, trajectory, sparse SfM point cloud and images of the City sequence.

Figure 4 :Figure 5 :

 45 Figure 4: Views of reconstructed surfaces of the city (real) sequence. Top: global and local views of the final surface (the sky is removed to help visualization). Bottom: local views of the surface during computation. The triangle normals are colored: the ground is white, the walls are red-green-blue, the sky is black.

Figure 6 :

 6 Figure 6: Images of the synthetic sequence, SfM results, the final and estimated surface.

  is maintained manifold by removing progressively tetrahedra from O t+1 until O t+1 ∩ D = / 0. Now we add P t+1 in T t without modifying neither O t+1 nor its border. If condition O t+1 ∩ D = / 0 can not be meet exactly, a minority of points in P t+1 are not added in T t to be sure that |δ O t+1 | is manifold. Last we trace rays to obtain the new freespace F t+1 , and grow O t+1 in F t+1 such that |δ O t+1 | is maintained manifold.

  .b). According to Sec. 3.3, we can use B = B r+l (c t+1 ).4.2 Step 2: O Shrinking Now we would like to find O t+1 such that O t+1 ⊆ O t \E (i.e. O t+1 ⊆ O t and O t+1 ∩E = / 0) and |δ O t+1 | is manifold. We initialize O t+1 = O t (Fig.

Acknowledgements

Thanks to CNRS, UBP, FEDER and Auvergne Région for funding. Thanks to the CRISTAL project for providing the ground truth VRML model used in the quantitative evaluation.