Adrien Bartoli
email: adrien.bartoli@inria.fr

Towards Gauge Invariant Bundle Adjustment: A Solution Based on Gauge Dependent Damping

Bundle ajustment is used to obtain accurate visual reconstructions by minimizing the reprojection error. The coordinate frame ambiguity, or more generality the gauge freedoms, has been dealt with in different manners. It has often been reported that standard bundle adjustment algorithms were not gauge invariant: two iterations within different gauges can lead to geometrically very different results. Surprisingly, most algorithms do not exploit gauge freedoms to improve performances. We consider this issue. We analyze theoretically the impact of the gauge on standard algorithms. We show that a sufficiently general damping matrix in Levenberg-Marquardt iteration can be used to implicitly reproduce a gauge transformation. We show that if the damping matrix is chosen such that the decrease in the reprojection error is maximized, then the iteration is gauge invariant. Experimental results on simulated and real data show that our gauge invariant bundle adjustment algorithm outperforms existing ones in terms of stability.

Introduction

The recovery of accurate 3D structure and camera motion from images is a major research challenge in photogrammetry [START_REF] Atkinson | Close Range Photogrammetry and Machine Vision[END_REF] and computer vision [START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF]. Bundle adjustment is a technique to refine a visual reconstruction to produce jointly optimal 3D structure and camera motion. Under certain assumptions on the noise on the observed features, bundle adjustment consists in minimizing the reprojection error, which is in general a non-linear procedure. Second-order non-linear least squares algorithms are usually employed, namely the Gauss-Newton and Levenberg-Marquardt methods. These methods iteratively improve sub-optimal parameter estimates by solving normal equations. Efficient solutions are possible thanks to the sparse block structure of the normal equations. The Levenberg-Marquardt method has proved the most successful due to the use of a trust region strategy, implemented via a damping of the normal equations, also called augmentation.

An inherent problem in bundle adjustment is the choice of the coordinate frame in which the reconstruction is expressed. This coordinate frame is called the gauge, or the datum in the photogrammetry community. A gauge is a subset of parameter vectors such that any two of them do not share the same underlying geometry. Parameter vectors corresponding to the same geometry are related by gauge transformations. It has been shown in [START_REF] Kanatani | Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy[END_REF] that these transformations have a group structure. The reprojection error is gauge invariant since it reflects the inherent merit of a reconstruction. These considerations hold whichever camera model is used and whichever calibration level is available. When ignored, gauge freedoms may induce numerical estimation problems, since they imply the rank deficiency of the normal equations.

It has been reported that most bundle adjustment algorithms [START_REF] Atkinson | Close Range Photogrammetry and Machine Vision[END_REF][START_REF] Dermanis | The photogrammetric inner constraints[END_REF][START_REF] Kanatani | Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy[END_REF][START_REF] Mclauchlan | Gauge invariance in projective 3D reconstruction[END_REF][START_REF] Mclauchlan | Gauge independence in optimization algorithms for 3D vision[END_REF][START_REF] Morris | Gauge fixing for accurate 3D estimation[END_REF][START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF] are not gauge invariant 1 . Carefully choosing the gauge can greatly improve the reliability. In practice, most algorithms converge to the same solution, but require different number of iterations. Surprisingly, gauge freedoms have rarely been exploited to improve on the reliability of algorithms.

Gauge freedoms imply that the normal equations to be solved at each iteration do not have a unique solution. Damping the normal equations consists in adding a symmetric positive definite matrix to the design matrix. This matrix is most of the time chosen as a diagonal one, which corresponds to an elliptical trust region, see e.g. [START_REF] Seber | Non-linear regression[END_REF] and §3 for more details.

The main contribution of this paper is an algorithm that maximizes the decrease in the reprojection error at each iteration, in a gauge invariant manner. In §4, we formally state the gauge invariance property of a bundle adjustment algorithm. In §5, we concentrate on Levenberg-Marquardt iterations and propose implicit gauges. In §6, we propose a gauge dependent damping matrix which (i) makes the underlying iteration gauge invariant, (ii) maximizes the decrease in the reprojection error up to first order and (iii) preserves the sparse block structure of the normal equations. Finally §7 reports some experimental results and §8 concludes. For simplicity of notation, we deal with projective reconstruction of points and cameras. Extension to other camera models, calibration levels, and types of features is straightforward, following [START_REF] Kanatani | Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy[END_REF][START_REF] Mclauchlan | Gauge independence in optimization algorithms for 3D vision[END_REF]. Note that an extended version of this paper is available.

Notation. We make no formal distinction between coordinate vectors and physical entities. Equality up to a non-null scale factor is denoted by ∼ and is sometimes written by introducing explicitly the scale factor: (x ∼ x ) ⇔ (∃α|x = αx ). Transposition and transposed inverse are denoted by T and -T . Vectors are typeset using bold fonts (q, Q), matrices using sans-serif fonts (P, T) and scalars in italics. Indices are used to indicate the size of a matrix or vector (P (3×4) , q (3×1) ) or to index a set of entities. The row-wise matrix vectorization is written vect.

Let P i , i = 1 . . . n denote the n reconstructed camera matrices and Q j , j = 1 . . . m denote the m reconstructed 3D points. Structure and motion parameters are contained in a (p × 1) vector X , partitionned into 12n motion parameters M and 4m structure parameters S as:

X T = (vect T (P 1 ) . . . vect T (P n ) M T Q T 1 . . . Q T m S T )
Changing the 15-degrees of freedom projective reconstruction basis and the individual scale factors of the camera matrices and point coordinate vectors leave the underlying geometry invariant. Hence there are g = 15 + n + m degrees of gauge freedom.

Previous work

A survey of gauge freedom handling strategies is [13, §9]. Methods to deal with this problem are either globally free, i.e. they left the gauge free to drift, or globally fixed, i.e. they enforce a global gauge. Globally free methods are based on selecting a solution of the normal equations, using a pseudo-inverse [START_REF] Atkinson | Close Range Photogrammetry and Machine Vision[END_REF], numerical damping [START_REF] Seber | Non-linear regression[END_REF][START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF] or local gauge constraints, for example, the standard photogrammetric inner constraints. Globally fixed methods enforce a prechosen gauge. Some methods use reference elements (trivial gauge), either features [START_REF] Atkinson | Close Range Photogrammetry and Machine Vision[END_REF][START_REF] Faugeras | What can be seen in three dimensions with an uncalibrated stereo rig?[END_REF] (object-centered gauge) or cameras [START_REF] Bartoli | On the non-linear optimization of projective motion using minimal parameters[END_REF][START_REF] Hartley | Euclidean reconstruction from uncalibrated views[END_REF][START_REF] Hartley | Projective reconstruction and invariants from multiple images[END_REF]. Another possibility for a globally fixed gauge is to use global gauge constraints [START_REF] Kanatani | Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy[END_REF][START_REF] Mclauchlan | Gauge invariance in projective 3D reconstruction[END_REF][START_REF] Mclauchlan | Gauge independence in optimization algorithms for 3D vision[END_REF], which raises the problem of enforcing these constraints during optimization. Besides the globally free photogrammetric inner constraints, all these methods are based on somehow arbitrary gauges which are not chosen to meet some efficiency criteria. A recent attempt in this direction is [START_REF] Morris | Gauge fixing for accurate 3D estimation[END_REF]. Based on gauge theory, the authors analyze, in the context of Euclidean reconstruction, which physical measure is the most likely to maximize the accuracy of the reconstruction while setting the scene scale.

Standard Algorithms

Bundle adjustment often consists in minimizing the reprojection error, given by the sum of squared differences between observed features q ij and predicted features q ij :

C(X ) = n i=1 m j=1 w ij d 2 (q ij , q ij ) = r T r, (1)
where w ij is one if point j is visible in view i and zero otherwise. The (r × 1) vector r is the residual error vector. Predicted features are given by q ij ∼ P i Q j . Given a sub-optimal solution X 0 , the parameter vector X is iteratively updated by X ← X + δ, where the increment δ is obtained as follows. Let J = ∂r ∂X denotes the (r × p) Jacobian matrix of the residual error vector r with respect to structure and motion parameters X , g = J T r be the (p × 1) gradient vector of C and let N = J T J be the Gauss-Newton approximation of the (p × p) Hessian matrix. The reprojection error is approximated by:

C(X + δ) ≈ C(X ) + g T δ + 1 2 δ T Nδ e . (2)
The minimum of this simple local quadratic model can be found by setting ∂e ∂δ = 0, which gives the Gauss-Newton iteration, through the following normal equations:

Nδ = -g. (3)
Note that N has a rank deficiency of g (the number of gauge freedoms). The Levenberg-Marquardt iteration is based on damping or augmenting the normal equations, as follows :

(N + W(λ) N )δ = -g, (4)
where λ > 0 is related to the trust region radius and W(λ) is some symmetric positive definite (p × p) weight matrix, called the damping matrix, often chosen as W(λ) = λI (p×p) , which corresponds to a spherical trust region. This is the original strategy proposed by Levenberg and Marquardt and recommended in [START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF]. Another commonly used solution, due to [START_REF] Seber | Non-linear regression[END_REF] and recommended in [START_REF] Hartley | Euclidean reconstruction from uncalibrated views[END_REF][START_REF] Press | Numerical Recipes in C -The Art of Scientific Computing[END_REF] is W(λ) = λD where D is a diagonal matrix containing the diagonal entries of N, which gives an elliptical trust region. The damping matrix must satisfy a normalization constraint so that the trust region radius is meaningful, e.g. I (p×p) = √ p. Note that the damping guarantees that N is full-rank. Parameter λ is tuned as follows: if parameters X + δ decrease the error, i.e. if C(X + δ) < C(X ), then the step is accepted and the value of λ is divided by some constant, often 10, else the step is rejected and λ is multiplied by the constant.

Influence of the Gauge

We formulate the gauge invariance property and examine under which conditions standard algorithms satisfy it.

Gauge Transformation and Invariance

Gauge transformations change the parameter vector without changing the underlying geometry, e.g. [START_REF] Kanatani | Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy[END_REF]. Let X T = (M T S T ) and X T = ( MT ŠT ) be two parameter vectors of the same reconstruction expressed within two different gauges G and Ǧ. Let T be the full-rank transformation relating the two underlying projective bases, defined such that Pi = γ i P i T and Qj = α j T -1 Q j , where γ i and α j are unknown non-zero scale factors. Entities written with a ˇare expressed within the gauge Ǧ. The structure and motion parameters transform as:

Š = diag(α 1 T -1 , . . . , α m T -1 m ) S M = diag(γ 1 T T , . . . , γ n T T 3n ) M,
since vect(γPT) = diag(γT T , γT T , γT T ) vect(P). We deduce that the parameter vectors are related by X = TX where the gauge transformation T is defined by:

T = diag(γ 1 T T , . . . , γ n T T 3n , α 1 T -1 , . . . , α m T -1 m ). (5)
Note that det(T) = 0 ⇒ det( T) = 0 since the γ i and α j are non-zero. In general, X and X are related by a unique gauge transformation.

The geometric equivalence of two reconstructions, denoted by ∼ =, is an egality 'up to the gauge' defined by:

(X ∼ = X ) ⇔ ∃ T | X = TX . (6)
with T of the form [START_REF] Hartley | Euclidean reconstruction from uncalibrated views[END_REF]. We can now formally state the gauge invariance property. A bundle adjustment algorithm is gauge invariant if it preserves geometric equivalence: (X ∼ = X ) ⇒ ((X + δ) ∼ = ( X + δ)). By writing explicitly this property using equation ( 6), expanding, and since X and X are related by a unique T, we obtain the following definition of gauge invariance:

∃ T | ( X = TX ) ⇒ ( δ = Tδ). (7)
Note that although derived in a different manner, this definition is similar to the one given in [START_REF] Mclauchlan | Gauge independence in optimization algorithms for 3D vision[END_REF].

Gauss-Newton

Consider a Gauss-Newton iteration. The increment δ within gauge G is given by solving the normal equations (3). Due to gauge freedoms, matrix N has a rank deficiency of g, and multiple solutions hold, corresponding to different gauges. Let G be a (p × g) full column-rank matrix defined by NG = 0, i.e. the columns of G span the nullspace of N. Denoting by † the Moore-Penrose pseudo-inverse, the possible solutions of the normal equations are parameterized by a (g × 1) vector v as δ = -(J T J)

† J T r + Gv, where we substituted N = J T J and g = J T r.

Since J = ∂r ∂X = ∂r ∂ X ∂ X ∂X = J T, we obtain: δ = -( T TJTJ T) † T TJ r + T -1 Ǧv. From equa- tion (7), gauge invariance holds if and only if ( T T Ň T) † = T -1 Ň † T -T and v = v
, which is verified if and only if T is orthonormal [START_REF] Mclauchlan | Gauge independence in optimization algorithms for 3D vision[END_REF], i.e. T -1 = T T . Hence, Gauss-Newtonbased bundle adjustment is not gauge invariant. The previously derived condition does not leave enough flexibility to be exploited in this direction.

Levenberg-Marquardt

Consider a Levenberg-Marquardt iteration. The increment δ within gauge G is given by solving the damped normal equations ( 4) N δ = -g. By substituting N = N + W(λ) = J T J + W(λ) and g = J T r in this equation and since det(N ) = 0, we obtain: δ = -(J T J + W(λ))

-1 J T r.

By substituting J = J T and expanding, we get:

δ = -( T T JT J T + W(λ)) -1 T T JT r = -T -1 ( Ň + T -T W(λ) T -1 ) -1 JT r.
From equation ( 7), gauge invariance holds if and only if:

W(λ) = T -T W(λ) T -1 . (8)
The usual choices for matrices W(λ) and W(λ), e.g. W(λ) = W(λ) = λI, do not fulfill the above-derived property 2 . Hence, standard implementations of Levenberg-Marquardt-based algorithms are not gauge invariant. Equation ( 8) can be verified if an appropriate gauge dependent choice is made for the damping matrices. This is the cornerstone for the gauge invariant method proposed in the next section.

Explicit and Implicit Gauges

Standard bundle adjustment algorithms are not gauge invariant. Hence, given a parameter estimate, there must exist an optimal local gauge within which the iteration reduces the reprojection error better than within the others. The problem of finding this gauge is dealt with in the next section. In this section, we examine the relationship between the damping matrix and gauge transformations.

A commonly used solution to globally enforce a gauge, e.g. [START_REF] Kanatani | Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy[END_REF][START_REF] Mclauchlan | Gauge invariance in projective 3D reconstruction[END_REF][START_REF] Mclauchlan | Gauge independence in optimization algorithms for 3D vision[END_REF], is what we call explicit gauge fixing. It consists in explicitly expressing the reconstruction within the desired gauge before each iteration.

We propose implicit gauge fixing. Consider equation ( 8). On the one hand, it tells us that the algorithm is not gauge invariant. On the other hand, it shows that there exists a close link between the gauge and the damping matrix. Hence, it can be used to perform bundle adjustment within gauge Ǧ, while choosing W(λ) such that it behaves exactly as if is was conducted within gauge G. For example if W(λ) = λI and W(λ) = λ T -T T -1 , then running the algorithm within gauges G and Ǧ is strictly equivalent. We confirm in our experiments that explicit and implicit gauges give exactly the same results. Let A = T T T, then W(λ) = λ T -T T -1 is given by:

W(λ) = λ diag(γ 2 1 A -1 , . . . , γ 2 n A -1 3n , α 2 1 A, . . . , α 2 m A m ). ( 9
)
Proposition 1 (Gauge invariance conditions) Levenberg-Marquardt iterations can be made gauge invariant if the following two conditions are verified: (i) the form of the damping matrix is at least as general as the form given by equation ( 9) and (ii) the damping matrix is uniquely determined by a gauge dependent criterion defined such that equation ( 8) is verified.

Gauge Dependent Damping

We propose a solution to choose a gauge dependent damping matrix W(λ) = λW which has the three following important properties. First, it makes the underlying iteration gauge invariant. Second, it maximizes the decrease in the reprojection error up to first order. Third, it preserves the sparse block structure of the normal equations.

In order to achieve this result, we begin by deriving an approximation of the reprojection error, depending upon W(λ). A quadratic approximation of C(X + δ) is given by equation ( 2). The incremental error e is given by e = g T δ + 1 2 δ T Nδ. The increment vector δ is given by solving the augmented normal equations (4) δ = -(N + W(λ))

-1 g.

Since det(W(λ)) = 0, we can formulate the following Taylor expansion

(N + W(λ)) -1 = ∞ n=0 (-1) n Z(n), where Z(n) = (W(λ) -1 N) n W(λ) -1 .

This expression is valid provided W(λ)

-1 N < 1, i.e. when N is small or λ is large.

This leads to δ = -∞ n=0 (-1) n Z(n)g. We substitute the above expression in equation [START_REF] Bartoli | On the non-linear optimization of projective motion using minimal parameters[END_REF]. The first term is rewritten as

g T δ = - ∞ n=0 (-1) n K(n), where K(n) = g T Z(n)g.
The second term becomes:

1 2 δ T Nδ = 1 2 ∞ n=0 (-1) n g T Z(n) N ∞ n=0 ((-1) n Z(n)g) .
After expansion and using the following property arising from the symmetry of W(λ): 1) n (n + 1)K(n + 1), which gives the following Taylor expansion for e:

Z(n)NZ(n ) = Z(n + n + 1), we obtain 1 2 δ T Nδ = 1 2 ∞ n=0 (-
e = ∞ n=0 (-1) n+1 1 + n 2 K(n). ( 10
)
Proposition 2 (Gauge invariant iterations) The damping matrix defined such that equation ( 10) is minimized satisfies proposition 1. Hence, the underlying iteration is gauge invariant.

Proof. We drop the parameter λ for this proof. Equation ( 10) can be rewritten as:

e = g T -W -1 + 3 2 W -1 NW -1 -. . . g.
Transfer the reconstruction from the gauge G to Ǧ by applying the gauge transformation T. By substituting g = T T ǧ and N = T T Ň T, we obtain:

e = ǧT T -W -1 + 3 2 W -1 T T Ň TW -1 -. . . T T ǧ,
and hence, W-1 = TW -1 T T , which verifies the gauge invariance equation ( 8) and concludes the proof.

Finding the damping matrix that minimizes e, equation ( 10), is a complicated problem. We keep only the first-order term: e ≈ -g T W(λ)

-1 g, which reduces the problem to finding a symmetric, positive definite, damping matrix λW W = arg max W, W 2 =p g T (λW) -1 g. The normalization condition W 2 = p is the same as the standard choice

W(λ) = λI ( I (p×p) 2 = p). Since W(λ) is a symmetric matrix, W(λ) -1 = 1 λ R T R,
where R is an upper-triangular matrix. With this parameterization, and since g T 1 λ R T Rg = 1 λ Rg 2 , we obtain:

R = arg max R, R T R 2 =p 1 λ Rg 2 . (11)
Hence, the coefficients of matrix R can be found by solving a simple linear least squares optimization problem.

From proposition 1 and imposing the fact that the damping should not spoil the sparse block structure of the normal equations, we are left with the following possible form:

R = diag(R M,1 , . . . , R M,n n , R S,1 , . . . , R S,m m ),
where the R M,i and the R S,j are respectively (12 × 12) and (4 × 4) upper-triangular matrices. Given that each of the p parameters gives one constraint on the damping matrix and since it must be uniquely determined (i.e. it must not have more than p parameters), we propose the following 10(n + 1)-parameter choice: R M,i = diag(M i , M i , M i ) and R S,j = S, where the M i and S are (4 × 4) upper-triangular matrices. This choice means that the parameter block of each camera i has its own 10 parameter damping block M T i M i , while a unique 10 parameter damping block S T S is defined for point parameters. This form is handy since each block can be found independently by solving a problem of the 1. Compute the damping matrix W(λ) -1 = λR T R where:

R = diag(M 1 , . . . , M n 3n , S, . . . , S m
).

The (4 × 4) upper-triangular matrices M i and S are formed by solving (see main text):

M i = arg max Mi, M T i Mi 2 =12 3 k=1 M i g k M,i 2 S = arg max S, S T S 2 =4 m j=1 Sg S,j 2 ,
where g S,j and g k M,i are (4×1) gradient vectors for the j-th point and for the k-th row of the i-th camera matrix respectively. 2. Perform one Levenberg-Marquardt iteration with the damping matrix W(λ). This iteration is gauge invariant due to the above choice. 3. Optional: Project the estimate onto a global gauge. Table 1. The gauge invariant iteration we propose. Periodical enforcement of global gauge constraints (step 3) such as renormalization of homogeneous coordinates is recommended. Gauge invariance holds in the sense that any global gauge can be enforced in step 3 without affecting the result. form [START_REF] Press | Numerical Recipes in C -The Art of Scientific Computing[END_REF]. More details are given in table 1. Below, we give details about solving the low-dimensional maximization problems of table 1. Let the problem to be solved be max x, x 2 =b Cx 2 . The singular vector v associated to the largest singular value of matrix C gives the solution for x as x = √ bv. Singular value decomposition of matrix C can be used to compute v.

Experimental Results

We compare our gauge invariant algorithm, denoted GAUGE INVARIANT, to some other ones. FREE directly optimizes the camera matrices and 3D points. Pseudo-inverse is used to solve the normal equations. HARTLEY [START_REF] Hartley | Projective reconstruction and invariants from multiple images[END_REF] consists in using a reference camera to partially eliminate the gauge and BARTOLI [START_REF] Bartoli | On the non-linear optimization of projective motion using minimal parameters[END_REF] consists in using two reference cameras to eliminate the gauge. FAUGERAS [START_REF] Faugeras | What can be seen in three dimensions with an uncalibrated stereo rig?[END_REF] is based on fixing the coordinates of 5 reference points to enforce the gauge.MCLAUCHLAN [START_REF] Mclauchlan | Gauge invariance in projective 3D reconstruction[END_REF] consists in enforcing a normalized basis after each iteration and in using first order gauge constraints in the normal equations and KANATANI [START_REF] Kanatani | Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy[END_REF] is similar to method FREE but periodically projects the estimate on a pre-chosen gauge to prevent it to drift too much. The initial solution, denoted by INIT, is computed by registering cam-eras in turn to a two-view reconstruction.

Simulated Data

We simulate m = 100 points lying in a cube with 1 meter side length, observed by n = 5 cameras with a focal length of 1000 pixels. The points are offset from a base plane lying inside the cube, with a mean offset denoted by d. Cameras are situated 10 meters away from the center of the cube. The baseline between consecutive cameras is 3 meters. All points are visible in all views. We add a centered Gaussian noise on true point positions with a 1 pixel variance. The results are averaged over 50 trials.

Figure 1 (left and middle) shows the results when the scene flatness, i.e. the mean offset d from the plane, is varied. We observe that the reprojection errors, i.e. the accuracy, are undistinguishable for all methods, except for methods FREE which converges to a different local minimum than the others, for weak geometry, i.e. when the mean offset from the plane d is small. Concerning the number of iterations, i.e. the computational cost, we observe that when the mean offset is large, i.e. when the geometry is strong, there is only slight differences between the different methods besides for method FREE, which takes clearly more iterations to converge. As the geometry becomes weaker, i.e. when the offset decreases, large discrepancies can be observed between the different methods. Method FAUGERAS based on reference points gives bad results, since using reference points to fix a projective basis can be very unstable. Methods HARTLEY and BARTOLI based on reference cameras and KANATANI based on periodical projection on a pre-chosen gauge give reliable results when the geometry is strong enough. Method MCLAUCHLAN performs better, while GAUGE INVARIANT is the method the less sensitive to the instability of the scene.

We tested other strong to weak scene configurations, based on varying the number of points and cameras, the baseline between consecutive cameras and the visibility (not shown here due to lack of space). We observed similar results as above: strong geometries reduce the discrepancies between the different methods. Similarly, we tested the influence of the initialization. As expected, discrepancies between the different methods reduce as the initialization gets more accurate.

Figure 1 (right) shows the reprojection error as a function of the iterations for the scene setting with d = 0.5. As expected, we observe that the error is faster to decrease as the number of iterations gets low.

Real Data

We compare the algorithms on various image streams. For two of them, the office sequence and the hotel sequence 3 sequence, they all gave 0.9594 pixel, besides method FREE which gave 1.5268 pixels. Concerning the number of iterations, the same observations as for simulated data can be made. The hotel sequence gives a weak geometry since an affine camera model is well-adapted to these data. Hence, a full perspective projection model is not well-constrained, which tends to increase the discrepancies between the different methods.

Algorithm

Conclusions

We proposed a bundle adjustment algorithm which maximizes the decrease in the error at each iteration, regardless the gauge within which the reconstruction is expressed. We derived this algorithm based on a careful study of how the gauge influence the iterations of standard algorithms, which are not gauge invariant. In particular, concerning Levenberg-Marquardt iterations, we showed that the standard diagonal damping matrices, defining elliptical trust regions, transform to a more complex shaped trust region when changing the gauge. Based on this, we proposed a gauge dependent damping matrix and a practical algorithm to compute it, that allows gauge invariant iterations while maximizing the decrease in the reprojection error. The sparse block structure of the normal equations is preserved.

We compared our algorithm to existing ones using simu-lated and real data. We observed that it is more reliable in the sense that it is less sensitive to weak geometry and weak initialization. Moreover, the error decreases more intensely throughout the iterations. Most algorithms converge to the same solution, but require different number of iterations.

Figure 1 .

 1 Figure 1. Reprojection error (left) and number of iterations (middle) when varying the scene unflatness from weak to strong geometry and reprojection error as a function of the iterations (right).

  , we show results, see table 2. For the office sequence, all methods give a reprojection error of 0.8392 pixel. For the hotel

		3												70
	frag replacements	2.5												60
	reprojection error (pixel)	0.5 1 1.5 2						INIT FREE FAUGERAS MCLAUCHLAN KANATANI BARTOLI HARTLEY		number of iterations	10 20 30 40 50
								GAUGE INVARIANT
		0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
							scene unflatness										scene unflatness

Table 2 .

 2 office seq. hotel seq. Number of iterations of the algorithms for the office and the hotel sequences.

	FREE	12	28
	HARTLEY	7	16
	BARTOLI	8	19
	FAUGERAS	10	22
	MCLAUCHLAN	6	15
	KANATANI	7	16
	GAUGE INVARIANT	7	13

We employ 'gauge invariant' in a manner different from[START_REF] Kanatani | Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy[END_REF][START_REF] Mclauchlan | Gauge invariance in projective 3D reconstruction[END_REF][START_REF] Mclauchlan | Gauge independence in optimization algorithms for 3D vision[END_REF]. These authors propose algorithms based on using a pre-defined global gauge. Hence, no matter the gauge within which the initial solution is expressed, their algorithms will give the same result. However, the behaviour of these algorithms do depend upon the pre-defined global gauge. In this sense, these algorithms are not gauge invariant.

It can be shown that equation[START_REF] Mclauchlan | Gauge invariance in projective 3D reconstruction[END_REF] with W(λ) = W(λ) = λI is verified as in the Gauss-Newton case, i.e. if T is orthonormal.

provided by the Modeling by Videotaping group in the Robotics Institute, Carnegie Mellon University.


Hyeong-Ohk Bae
email: hobae@ajou.ac.kr

Lorenzo Brandolese
email: brandolese@math.univ-lyon1.fr

On the effect of external forces on incompressible fluid motions at large distances

Keywords: Navier-Stokes, asymptotic profiles, asymptotic behavior, far-field, spatial infinity, weak solutions, strong solutions, incompressible viscous flows, decay estimates. 2000 Mathematics subject classification 76D05, 35Q30

Introduction

Let d ≥ 2 be an integer. We consider the Navier-Stokes equations in R d , for an incompressible fluid submitted to an external force f :

R d × R + → R d :      ∂ t u + u • ∇u = ∆u -∇p + f divu = 0 u(x, 0) = a(x) (1)
We want to investigate the effect of the external force on the large time and spatial asymptotics of the solution. In particular, we will show that even if f is small and well localized (say, compactly supported in space-time) but with non-zero mean, then the velocity of the fluid particles at all times t > 0 and in all point x outside of balls B(0, R(t)) with large radii, is consirerably faster than in the case of the free Navier-Stokes equations.

This sharply contrasts with the asymptotic properties of the solution of other semilinear parabolic equations, where localized external forces do not affect the behavior of the solution as |x| → ∞: only the large time behavior is influenced by the force.

The main results of the paper will be stated in Section 3: they include new sharp pointwise estimates of the form

ct|x| -d ≤ |u(x, t)| ≤ c ′ t|x| -d , (2)
valid for 0 < t < t 0 small enough and |x| ≥ R(t). Such pointwise estimates should be compared to those in the case of the Navier-Stokes equations without forcing studied in [START_REF] Brandolese | New Asymptotic Profiles of nonstationnary solutions of the Navier-Stokes system[END_REF], where the behavior of u is like t|x| -d-1 . The constant c in (2) essentially depends on the integral of f . Estimates (2) can be applied, e.g. to the case of compactly supported initial data. In such case, they describe how fast fluid particles start their motion in the far field, at the beginning of the evolution. An even more precise description of the motion of flows with localized data at large distances will be given by the following asymptotic profile, valid for all t > 0 such that the strong solution u is defined and for all |x| ≥ R(t):

u(x, t) ≃ K(x) t 0 f (y, s) dy ds, (3)
where K(x) is the matrix of the second order derivatives of the fundamental solution of the Laplacian in R d . As such, K j,k (x) is a homogeneneous function of degree -d. See Theorem 1, Section 3, for a more precise statement. As a consequence of our asymptotic profiles one can also recover some known bounds on the large time behavior of L p -norms for strong solutions with very short proofs. See, e.g., [START_REF] Amrouche | Schonbek Pointwise Decay of Solutions and of higher derivatives to Navier-Stokes Equations[END_REF][START_REF] Bae | Upper and lower bounds of temporal and spatial decays for the Navier-Stokes equations[END_REF][START_REF] Choe | Weighted estimates of the asymptotic profiles of the Navier-Stokes flow in R n[END_REF][START_REF] He | On the decay properties of Solutions to the nonstationary Navier-Stokes Equations in R 3[END_REF][START_REF] Kukavica | Weighted L p decay for solutions of the Navier-Stokes equations[END_REF][START_REF] Oliver | Remark on the rate of decay of higher order derivatives for solutions to the Navier-stokes equations in R n[END_REF][START_REF] Schonbek | L 2 decay for Weak Solutions of the Navier-Stokes Equations[END_REF][START_REF] Schonbek | Lower Bounds of Rates of Decay for Solutions to the Navier-Stokes Equations[END_REF][START_REF] Skaák | Asymptotic decay of higher-order norms of solutions to the Navier-Stokes equations in R 3[END_REF][START_REF] Qi | Global solutions of the Navier-Stokes equations with large L 2 norms in a new function space[END_REF][START_REF] Zhou | A remark on the decay of solutions to the 3-D Navier-Stokes equations[END_REF] for a small sample of recent works on this topic and related developements. Our approach is different from that of these papers since it consists in deducing information on spatial norms about the solution from information on their pointwise behavior. As such, it is not so well suited for the study of weak solutions as it is for that of strong, small solutions. However, its advantage is that it allows to get sharp estimates, both from above and below, also in the case of L p -norms with weight like, e.g., (1 + |x|) α , for a wider range of the parameters p and α. Namely, we will prove that, under localization assumptions on the datum and the force, and provided

∞ 0 f (y, s) dy ds = 0 R d , then, for 1 < p ≤ ∞, α ≥ 0 and t → ∞, (1 + |x|) α u(t) p ≃ ct -1 2 (d-α-d/p) , if α + d/p < d.
The above restriction on the parameters p and α are optimal, as we will show that, for all t > 0, and 1 ≤ p < ∞,

(1 + |x|) α u(t) p = +∞ if α + d/p ≥ d.
See Theorem 2, Section 3, for more precise statements.

Preliminary material

The assumptions on the external force will be the following:

|f (x, t)| ≤ ǫ (1 + |x|) -d-2 ∧ (1 + t) -(d+2)/2 (4) f L 1 (R d ×R + ) ≤ ǫ (5)
for some (small) ǫ > 0 and a.e. x ∈ R d and t ≥ 0. Note that assumption (5) can be interpretated as a logarithmic improvement on the decay estimate [START_REF] Brandolese | Asymptotic behavior of the energy and pointwise estimates for solutions to the Navier-Stokes equations[END_REF]. Of course, an additional potential force ∇Φ could also be added, affecting in this way the pressure of the fluid, but not its velocity field. On the other hand, because of condition [START_REF] Brandolese | New Asymptotic Profiles of nonstationnary solutions of the Navier-Stokes system[END_REF], and due to the unboundedness of singular integrals (the Riesz transforms) in L 1 and, we are not allowed to apply the Helmoltz decompsition to f . Therefore, we will not restrict ourselves to divergence-free external forces, as it is sometimes the case in the literature.

The above assumptions on f , as well as the smallness assumption on the datum below, look stringent: we simply put assumptions that allow us to provide the shortest possible and self-contained construction of strong solutions with some decay. Indeed, our main goal will be to show that even in the case of very nicely behaved external forces (with non-zero mean) and data, the solution will be badly behaved at infinity. More precisely, upper bounds on f and a = u(0), no matter how good, will lead to lower bounds on u(t).

But for the time being, we concentrate on the simpler problem of upper bounds on u, and start by establishing the following simple result.

Proposition 1 Assume that f satisfies (4)- [START_REF] Brandolese | New Asymptotic Profiles of nonstationnary solutions of the Navier-Stokes system[END_REF]. Let also a ∈ L 1 (R d ) be a divergence-free vector field such that, a 1 < ǫ and |a(x)| ≤ ǫ(1 + |x|) -d for a.e. x ∈ R d . If ǫ > 0 is small enough, then there exists a unique strong solution u of the Navier-Stokes equation (NS) satisfying, for some constant C > 0, the pointwise decay estimates

|u(x, t)| ≤ Cǫ (1 + |x|) -d ∧ (1 + t) -d/2 , (6)
and such that u(0) = a (in the sense of the a.e. and distributional convergence as t → 0).

In particular, for α ≥ 0, 1 < p ≤ ∞ , we have

(1 + |x|) α u(t) p ≤ C(1 + t) -1 2 (d-α-d p ) , when α + d/p < d. ( 7
)
The above estimate remains true in the limit case (α, p) = (d, ∞).

Remark 1

The Proposition above can be viewed as the limit case of a previous result by Takahashi [START_REF] Takahashi | A weighted equation approch to decay rate estimates for the Navier-Stokes equations[END_REF], where pointwise decay estimates of the form |u(x,

t)| ≤ C (1 + |x|) -γ ∧ (1 + t) -γ/2
] had been obtained for 0 ≤ γ < d, with a different method and under different assumptions.

There are, on the other hand, many other methods for proving sharp upper bounds of the form

(1 + |x|) α u(t) p ≤ Ct -1 2 (d-α-d/p) , .
For example, it would be possible to adapt the arguments of [START_REF] Choe | Weighted estimates of the asymptotic profiles of the Navier-Stokes flow in R n[END_REF] or [START_REF] Kukavica | Weighted L p decay for solutions of the Navier-Stokes equations[END_REF] to the case of non-zero external forces. See also [START_REF] Amrouche | Schonbek Pointwise Decay of Solutions and of higher derivatives to Navier-Stokes Equations[END_REF], [14] [19] for other different approaches for getting decay estimates. However, the optimal range of the parameters α and p for the validity of such estimate is not discussed in the previous papers. Our condition α + d/p < d is more general, and turns out to be optimal whenever f has non-zero mean, as we will show in Theorem 2.

We denote by P be the usual Leray projector onto the solenoidal vector fields and with e t∆ the heat semigroup. The solution of Proposition 1 is obtained by a straightforward fixed point argument by solving the equivalent integral formulation:

     u(t) = e t∆ a - t 0 e (t-s)∆ P∇ • (u ⊗ u)(s) ds + t 0 e (t-s)∆ Pf (s) ds diva = 0. (8)
Let K(x, t) be the kernel of e t∆ P. It is well known, and easy to check that, for all t > 0, K(•, t) belongs to C ∞ (R d ) and satisfies the scaling relation K(x, t) = t -d/2 K(x/ √ t, 1), and the decay estimates

|K(x, t)| ≤ C |x| -d ∧ t -d/2 , ( 9
) |∇K(x, t)| ≤ C |x| -d-1 ∧ t -(d+1)/2 . ( 10
)
The following Lemma is an immediate consequence of the above decay estimates.

Lemma 1 For all (x, t) such that |x| ≥ e √ t and t > 0, the kernel K(x, t) satisfies:

t 0 |y|≤|x| |K(y, s)| dy ds ≤ Ct log(|x|/ √ t). (11)
The proof of Proposition 1 relies on the following lemma.

Lemma 2 Let f be such that (4)-( 5) holds. Let L be the linear operator defined by

L(f ) = t 0 e (t-s)∆ Pf (s) ds. ( 12
)
Then there is a constant C > 0 such that

|L(f )(x, t)| ≤ Cǫ (1 + |x|) -d ∧ (1 + t) -d/2 . ( 13
)
Proof. We give a proof based on Lorentz spaces L p,q , see [START_REF] Bergh | Interpolation Spaces, an Introduction[END_REF][START_REF] Lemarié-Rieusset | Recent developements in the Navier-Stokes problem[END_REF] for their definition and basic properties, since Lorentz spaces will play a role also in the sequel. For simplicity, we will drop the constant ǫ throughout the proof. By inequality (4) we readily get

f (t) L p,∞ ≤ C(1 + t) -1 2 (d+2-d p ) , 1 < p < ∞.
Now, interpolating the L p,q -space between L p 1 ,∞ and L ∞ , with 1 < p 1 < p and 1 ≤ q ≤ ∞, we deduce

f (t) L p,q ≤ C(1 + t) -1 2 (d+2-d p ) , 1 < p < ∞, 1 ≤ q ≤ ∞. (14)
In particular, choosing p = q = d,

f (s) L d ≤ C(1 + s) -(d+1)/2 .
On the other hand, it follows from (9) (or by well known L p -estimates on the kernel K),

that K(•, t -s) d/(d-1) ≤ C(t -s) -1/2 . By Young inequality we get K(t -s) * f (s) ∞ ≤ C(t -s) -1/2 (1 + s) -(d+1)/2 .
Another obvious estimate estimate is

K(t -s) * f (s) ∞ ≤ C(t -s) -d/2 f (s) 1 .
Now splitting the intergral defining L at s = t/2 we get

L(f )(t) ∞ ≤ C(1 + t) -d/2 .
In particular, for |x| ≤ e √ t , we obtain |L(f )(x, t)| ≤ C(1 + |x|) -d and now we can limit ourselves to the case |x| ≥ e √ t and |x| ≥ 1. We write

L(f ) = L 1 (f ) + L 2 (f ), where L 1 (f ) = t 0 |x-y|≤|x|/2 . . . and L 2 (f ) = t 0 |x-
y|≥|x|/2 . . . . For treating L 1 (f ), we use assumption (4) and obtain, after changing the variables z = x -y and τ = t -s, by Lemma 1,

|L 1 (f )(x, t)| ≤ Ct|x| -d-2 log(|x|/ √ t) ≤ C|x| -d ,
where we used that |x| ≥ e √ t. Using estimate [START_REF] Kato | Strong L p -Solutions of the Navier-Stokes Equations in R m , with applications to weak solutions[END_REF] and assumption (5

) we get |L 2 (f )(x, t)| ≤ C|x| -d .
Direct estimates on the heat kernel show that for a ∈ L 1 (R n ) satisfying a 1 ≤ ǫ and |a(x)| ≤ ǫ(1 + |x|) -d , we have, for some constant C > 0 independent on x and t ≥ 0,

|e t∆ a(x)| ≤ Cǫ (1 + |x|) -d ∧ (1 + t) -d/2 . ( 15
)
Owing to estimate [START_REF] Miyakawa | On space time decay properties of nonstationary incompressible Navier-Stokes flows in R n[END_REF], the fixed point argument argument used in Miyakawa [START_REF] Miyakawa | On space time decay properties of nonstationary incompressible Navier-Stokes flows in R n[END_REF], where the solution was constructed in the special case f ≡ 0, goes through in our situation. More precisely, the approximate solutions u k , constructed in the usual manner (see, e.g. [START_REF] Kato | Strong L p -Solutions of the Navier-Stokes Equations in R m , with applications to weak solutions[END_REF][START_REF] Miyakawa | On space time decay properties of nonstationary incompressible Navier-Stokes flows in R n[END_REF]) converge in the Banach space of measurable functions of the (x, t)-variables bounded by the right-hand side of [START_REF] Schonbek | L 2 decay for Weak Solutions of the Navier-Stokes Equations[END_REF]. The existence and the unicity of a solution u satisfying the pointwise estimates ( 6) follows.

It is straightforward to see that any function satisfying estimates (6) must also verify the weighted-L p bounds [START_REF] Choe | Weighted estimates of the asymptotic profiles of the Navier-Stokes flow in R n[END_REF]. This can be proved first by bounding the L p,∞ -norms of (1+|x|) α u(t) and then interpolating the L p -space between L p 1 ,∞ and L ∞ with 1 < p 1 < p. This is exactly the same argument that we applied to deduce estimate [START_REF] Oliver | Remark on the rate of decay of higher order derivatives for solutions to the Navier-stokes equations in R n[END_REF].

This completes the proof of Proposition 1.

Denote by Γ the Euler Gamma function and by E d the fundamental solution of -∆ in R d . The following Lemma (proved in [START_REF] Brandolese | Fine properties of self-similar solutions of the Navier-Stokes equations[END_REF]) will be useful:

Lemma 3 Let K = (K j,k ), where K j,k (x) is the homogeneous function of degree -d K j,k (x) = ∂ 2 x j ,x k E d (x) = Γ(d/2) 2π d/2 • -δ j,k |x| 2 + dx j x k |x| d+2 . ( 16
)
Then the following decomposition holds:

K(x, t) = K(x) + |x| -d Ψ x/ √ t , ( 17
)
where Ψ is a function defined on R d , smooth outside the origin and such that, for all α ∈ N d , and x = 0, |∂ α Ψ(x)| ≤ Ce -c|x| 2 , where C and c are positive constant, depending on |α| but not on x.

Before going further, let us put an additional assumption on the external force, namely

|x|f (•, t) 1 ≤ C(1 + t) -1/2 , for some C > 0 and all x ∈ R d , t ≥ 0. (18)
Just like condition (5), we can view assumption [START_REF] Takahashi | A weighted equation approch to decay rate estimates for the Navier-Stokes equations[END_REF] as another logarithmic improvement on the decay estimate (4). However, here the smallness of the constant C > 0 is unessential.

Next Lemma provides an explicit asymptotic expansion for L(f )(x, t), as |x| → ∞. The function of the (x, t)-variables R f defined by relation [START_REF] Qi | Global solutions of the Navier-Stokes equations with large L 2 norms in a new function space[END_REF] below can be viewed as a remainder term, that can be neglected at large distances, when |x| > > √ t .

Lemma 4 Let f such that assumptions (4)-( 5) and (18) hold. For x = 0, and t > 0, define R f (x, t), through the relation

L(f )(x, t) = K(x) t 0 f (y, s) dy ds + R f (x, t). ( 19
)
Then R f (x, t) satisfies, for some constant C > 0,

∀ |x| ≥ e √ t : |R f (x, t)| ≤ C|x| -d-1 √ t. ( 20
)
Remark 2 The product on the right hand side of ( 19) is the usual product between the d× d matrix K(x) and the d-vector field t 0 f k (y, s) dy ds. The j -component is thus given by d k=1 K j,k (x) t 0 f k (y, s) dy ds.

Proof. We introduce a function φ(x, t) through the relation

f (x, s) = f (y, t) dy g(x) + φ(x, t), (21)
where g t (x) = (4πt) -d/2 exp(-|x| 2 /(4t)) is the gaussian and g = g 1 . Identity ( 21) is inspired by a paper by Schonbek, [START_REF] Schonbek | Lower Bounds of Rates of Decay for Solutions to the Navier-Stokes Equations[END_REF]. See also [START_REF] Brandolese | New Asymptotic Profiles of nonstationnary solutions of the Navier-Stokes system[END_REF], where a similar idea was used to write an asymptotic expansion of the nonlinear term. From e t∆ g = g t+1 and from the fact that P commutes with the heat kernel we get (applying also Lemma 3) the decomposition

L(f )(x, t) = K(x) t 0 f (y, s) dy ds + R 1 (x, t) + R 2 (x, t),
where

R 1 (x, t) = |x| -d t 0 Ψ(x/ √ t + 1 -s) f (y, s) dy ds and R 2 (x, t) = t 0 K(x -y, t -s)φ(y, s) dy ds.
Using that |Ψ(x)| ≤ C|x| -1 (this is a consequence of Lemma 3) we obtain

|R 1 (x, t)| ≤ C|x| -d-1 t 0 √ t + 1 -s f (s) 1 ds ≤ C √ t|x| -d-1 .
For treating R 2 we use that φ(y, s) dy = 0 and apply the Taylor formula. We thus obtain the decomposition

R 2 = R 2,1 + • • • + R 2,4 , with R 2,1 = - t 0 |y|≤|x|/2 1 0 ∇ x K(x -θy, t -s) dθ • yφ(y, s) dy ds, next R 2,2 + R 2,3 = t 0 |x-y|≤|x|/2 + t 0 |y|≥|x|/2, |x-y|≥|x|/2
K(x -y, t -s)φ(y, s) dy ds, and

R 2,4 = - t 0 K(x, t -s) |y|≥|x|/2
φ(y, s) dy ds.

The term R 2,1 is bounded using estimate [START_REF] Lemarié-Rieusset | Recent developements in the Navier-Stokes problem[END_REF] and the estimate (deduced from ( 18))

|x|φ(s) 1 ≤ C(1 + s) -1/2 . This yields |R 2,1 (x, t)| ≤ C|x| -d-1 √ t.
The other three terms can be treated observing that

|φ(y, s)| ≤ C(1 + |y|) -d-2 .
Then applying Lemma 1 we get,

for |x| ≥ e √ t, |R 2,2 (x, t)| ≤ Ct|x| -d-2 log(|x|/ √ t) ≤ C|x| -d-1 √ t.
For the R 2,3 term, we can observe that the integrand is bounded by

C|y| -2d-2 . Therefore, |R 2,3 (x, t)| ≤ Ct|x| -d-2 ≤ C|x| -d-1 √ t.
The term R 2,4 can be estimated in the same way.

Main results

From the previous Lemma we now deduce the following result: it completes to the case f = 0 the asymptotic profile constructed in [START_REF] Brandolese | New Asymptotic Profiles of nonstationnary solutions of the Navier-Stokes system[END_REF].

Theorem 1 Let f , a and u be as in Proposition 1. We assume that f satisfies also condition [START_REF] Takahashi | A weighted equation approch to decay rate estimates for the Navier-Stokes equations[END_REF]. Then (the notation is explained in Remark 2)

u(x, t) = e t∆ a(x) + K(x) t 0 f (y, s) dy ds + R(x, t), ( 22
)
for some function R satisfying

∀ |x| ≥ e √ t : |R(x, t)| ≤ C|x| -d-1 √ t, ( 23
)
where C > 0 is a constant independent on x and t.

Remark 3 This theorem essentially states that if t > 0 is fixed and

t 0 f (y, s) dy ds = 0 R d , then u(x, t) ≃ e t∆ a(x) + K(x) t 0 f (y, s) dy ds as |x| → ∞.
Proof. Owing to our previous Lemma, the only thing that remains to do is to write u through Duhamel formula [START_REF] He | On the decay properties of Solutions to the nonstationary Navier-Stokes Equations in R 3[END_REF] and show that

B(u, u) ≡ t 0 e (t-s)∆ P∇ • (u ⊗ u)(s) ds can be bounded by C|x| -d-1 √ t.
In fact, a stronger estimate will be proved. Note that the convolution kernel F (x, t) of e t∆ P∇ satisfies to the same decay estimates of ∇K (see [START_REF] Lemarié-Rieusset | Recent developements in the Navier-Stokes problem[END_REF]). Therefore, after writing

B(u, u)(x, t) = t 0 F (x -y, t -s)(u ⊗ u)(y, s) dy ds,
then splitting the spatial integral into |y| ≤ |x|/2 and |y| ≥ |x|/2, and using that

F (•, t - s) 1 = c(t -s) -1/2 , we get |B(u, u)|(x, t) ≤ C|x| -d-1 t 0 u 2 2 ds + (1 + |x|) -2d √ t.
When d ≥ 3, by estimate ( 7) with α = 0 and p = 2, we get 2 2 ds ≤ C log(1 + t). (such facts on the L 2 -norm on u also follow from Schonbek's results [START_REF] Schonbek | L 2 decay for Weak Solutions of the Navier-Stokes Equations[END_REF][START_REF] Schonbek | Lower Bounds of Rates of Decay for Solutions to the Navier-Stokes Equations[END_REF]). This yields in particular the required bound

t 0 u(s) 2 2 ds ≤ C(1 ∧ t), which is enough to conclude. When d = 2, we have only t 0 u(s)
|B(u, u)|(x, t) ≤ C √ t|x| -d-1 .
Remark 4 The asymptotic profile (22) leads us to study the homogeneous vector fields of the form:

m(x) = K(x) c, c = (c 1 , . . . , c d ).
Notice that m has a zero in R d \{0} if and only if c = 0. Indeed, we can limit ourselves to the points ω ∈ S d-1 and recalling [START_REF] Schonbek | Lower Bounds of Rates of Decay for Solutions to the Navier-Stokes Equations[END_REF] we see that m vanishes at the point ω if and only if

dω j ω • c = c j , for j = 1, . . . , d.
Multiplying scalarly with ω we get (d -1)ω • c = 0 for all ω ∈ S d-1 . As d ≥ 2, we obtain c = 0. Applying this observation to our situation we get

inf x =0 |x| d K(x) t 0
f (y, s) dy ds > 0 for all t > 0 such that t 0 f (y, s) dy ds = 0.

A first interesting consequence of Theorem 1 is the following:

Corollary 1 1. Let a, f and u as in Theorem 1. We additionally assume that a = u(0) satisfies a(x) = o(|x| -d ) as |x| → ∞. Let t > 0 be fixed, such that t 0 f (y, s) dy ds = 0 R d . Then for some locally bounded function R(t) > 0 (R(t) can blow up for t → ∞ or t → 0) and all |x| ≥ R(t), we have

c t |x| -d ≤ |u(x, t)| ≤ c ′ t |x| -d (24)
for some constant c t , c ′ t > 0 independent on x 2. (Short time behavior of the flow)

If in addition f ∈ C(R + , L 1 (R d )) and f 0 = f (•, 0
) is such that f 0 (y) dy = 0 R d , then the behavior of u at the beginning of the evolution can be described in the following way: there exists a time t 0 > 0 such that for all 0 < t < t 0 and all |x| ≥ R(t)

ct|x| -d ≤ |u(x, t)| ≤ c ′ t|x| -d ,
for some constants c, c ′ > 0 independent on x and t.

Proof. Let

m t (x) ≡ K(x) t 0
f (y, s) dy ds.

We apply the asymptotic expansion ( 22) for u and study each term in this expression. By the non-zero mean condition on f and Remark 4, 3| m t (x)| can be bounded from above and from below as in (24), for some 0 < c t < c ′ t . Estimate (23) shows that |R(x, t)| ≤ c t |x| -d provided |x| ≥ R(t) and R(t) > 0 is taken sufficiently large.

On the other hand, the bound In conclusion the first and third terms in expansion ( 22) can be absorbed by m t (x). The first conclusion of the corollary then follows. The second conclusion is now immediate, because for t > 0 small enough and |x| ≥ R(t) large enough, one can find two constants α, β > 0 such that

|g t (x -y)| ≤ C|x -y| -d-
αt|x| -d ≤ | m t (x)| ≤ βt|x| -d
(the condition f 0 (y) dy = 0 R d is needed only for the lower bound).

Remark 5

The result for the free Navier-Stokes equation (f ≡ 0) is different (see [START_REF] Brandolese | New Asymptotic Profiles of nonstationnary solutions of the Navier-Stokes system[END_REF]): instead of (24) one in general obtains estimates of the form

c t |x| -d-1 ≤ |u(x, t)| ≤ c ′ t |x| -d-1
. Therefore, an external force, even if small and compactly supported in space-time, has the effect of increasing the velocity of the fluid particles at all points at large distances.

We obtain, as another consequence of Theorem 1, the following result for the large time behavior.

Theorem 2 Let a, f and u as in Theorem 1. We additionally assume that a = u(0) satisfies a

(x) = o(|x| -d ) as |x| → ∞ and that ∞ 0 f (y, s) dy ds = 0 R d .
Then, for α ≥ 0, 1 < p < ∞, for some t 0 > 0, some constants c, c ′ > 0 and and all t > t 0 we have,

ct -1 2 (d-α-d p ) ≤ (1 + |x|) α u(t) p ≤ c ′ t -1 2 (d-α-d p ) , when α + d/p < d. ( 26
)
The above inequalities remain true in the limit case (α, p) = (d, ∞).

On the other hand,

(1 + |x|) α u(t) p = +∞ when α + d/p ≥ d, (27)
and the above equality remains true in the limit cases p = 1 or (α > d and p = ∞).

Proof. The upper bound in (26) has been already proved in Proposition 1. The proof makes use of an argument used before in [START_REF] Brandolese | New Asymptotic Profiles of nonstationnary solutions of the Navier-Stokes system[END_REF]. By our assumption on f and Remark 4, the homogeneous function f (y, s) dy ds. We get for a sufficiently large M > 0, all t > M and all |x| ≥ M √ t,

m(x) ≡ K(x) ∞ 0 f (
|u(x, t) -e t∆ a| ≥ | m(x)| -C|x| -d-1 √ t -c 0 2 |x| -d ≥ c 0 3 |x| -d .
On for some c > 0 and all t > M . In the same way, inequality (28) also implies the lower bound (26) in the limit case p = ∞ and conclusion (27). The upper bounds for u obtained in Proposition 1 complete the proof.

The lower bounds obtained in (26) are invariant under the time translation t → t + t * , but the hypothesis made on f , ∞ 0 f (y, s) dy ds = 0 R d , is not. The explanation is that the conditions on the datum, a(x) = o(|x| -d ) as |x| → ∞, and a ∈ L 1 (R d ), in general, are not conserved at later times. In fact, according to ( 22 We did not treat the case of external forces with vanishing integrals. In that case, the principal term in the asymptotic expansion [START_REF] Qi | Global solutions of the Navier-Stokes equations with large L 2 norms in a new function space[END_REF] of u as |x| → ∞ disappears. A similar method, however, can be used to write the next term in the asymptotics which equals ∇K(x) : t 0 y ⊗ f (y, s) dy ds

(more explicitly, h,k ∂ h K j,k (x)
t 0 y h f k (y, s) dy ds, for the j -component, j = 1, . . . , d). For this we need to put more stringent assumptions on the decay of f : the spatial decay must be increased of a factor (1 + |x|) -1 and the time decay by a factor (1 + t) -1/2 . Assumptions ( 5)-( 18) should also be sharpened accordingly. For well localized initial data a(x), then one would deduce, bounds of the form

c t |x| -d-1 ≤ |u(x, t)| ≤ c ′ t |x| -d-1
and

(1 + |x|) α u(t) p ≃ ct -1 2 (d+1-α-d/p) , for α + d/p < d + 1.

However, as for the free Navier-Stokes equations (see [START_REF] Brandolese | New Asymptotic Profiles of nonstationnary solutions of the Navier-Stokes system[END_REF]), suitable additional non-symmetry conditions on the flow and on the matrix t 0 y ⊗ f (y, s) dy ds should be added for the validity of the lower bounds.

1 √ t and g t 1 = 1 2 √

 1112 for the heat kernel lead to |e t∆ a(x)| ≤ C |y|≤|x|/t|x -y| -d-1 |a(y)| dy + |y|≥|x|/2 g t (x -y)|a(y)| dy ≤ C|x| -d-1 √ t a 1 + ess sup |y|≥|x|/2 |a(y)| (25) and since |a(y)| = o(|y| -d ) as |y| → ∞, this expression is also bounded by c t |x| -d for |x| ≥ R(t) and a sufficiently large R(t).

  y, s) dy ds does not vanish. Therefore, for some c 0 > 0, we obtain | m(x)| ≥ c 0 |x| -d . Let us apply the profile (22), writing the second term on the right hand side as m(x)-K(x) ∞ t

  ), one has |u(x, t * )| = o(|x| -d ) as |x| → ∞ if and only if t * 0 f (y, s) dy ds = 0 R d .

  the other hand the computation (25) guarantees that we can bound, for |x| ≥ M √ t and t > M , |e t∆ a(x)| ≤ c 0 12 |x| -d . We thus get Multiplying this inequality by the weight (1 + |x|) αp , then integrating with respect to x on the set |x| ≥ M √ t , we immediately deduce (1 + |x|) α u(x, t)

	|u(x, t)| ≥ c 0 4 |x| -d ,	for all |x| ≥ M	√	t, t > M .	(28)
	Let 1 < p ≤ ∞. p	≥ ct -1 2 (d-α-d p )

Acknowledgments

The preparation of this paper was supported by EGIDE, through the program Huber-Curien "Star" N. 16560RK


Adrien Bartoli

The Geometry of Dynamic Scenes -On Coplanar and Convergent Linear Motions Embedded in 3D Static Scenes

In this paper, we consider structure and motion recovery for scenes consisting of static and dynamic features. More particularly, we consider a single moving uncalibrated camera observing a scene consisting of points moving along straight lines converging to a unique point and lying on a motion plane. This scenario may describe a roadway observed by a moving camera whose motion is unknown.

We show that there exist matching tensors similar to fundamental matrices. We derive the link between dynamic and static structure and motion and show how the equation of the motion plane (or equivalently the plane homographies it induces between images) may be recovered from dynamic features only.

Experimental results on real images are provided, in particular on a 60frames video sequence.

Introduction

Most existing works on the geometry of multiple images rely on the assumption that the observed scene is rigid. The rigidity constraint allows to derive matching relations among two or more images, represented by e.g. the fundamental matrix or trifocal tensors. These matching tensors encapsulate the motion and the intrinsic parameters of the cameras which took the underlying images, and thus all the geometric information needed to perform 3D reconstruction. Matching tensors for rigid scenes can also be employed for scenes composed of multiple, independently moving objects [START_REF] Costeira | A multi-body factorization method for motion analysis[END_REF][START_REF] Fitzgibbon | Multibody structure and motion: 3-D reconstruction of independently moving objects[END_REF][START_REF] Wolf | Two-body segmentation from two perspective views[END_REF], which requires however that enough features be extracted for each object, making segmentation, at least implicitly, possible.

On the other hand, there is a growing body of literature [START_REF] Avidan | Trajectory triangulation: 3D reconstruction of moving points from a monocular image sequence[END_REF][START_REF] Han | Reconstruction of a scene with multiple linearly moving objects[END_REF][START_REF] Han | Multiple motion scene reconstruction from uncalibrated views[END_REF][START_REF] Shashua | Homography tensors: On algebraic entities that represent three views of static or moving planar points[END_REF][START_REF] Sturm | Structure and motion for dynamic scenes -the case of points moving in planes[END_REF][START_REF] Wolf | On projection matrices P k → P 2 , k = 3, • • • , 6, and their applications in computer vision[END_REF][START_REF] Wolf | Join tensors: On 3D-to-3D alignment of dynamic sets[END_REF] dealing with the case of independently moving features, often termed as dynamic features. The goal of these works is to provide algorithms for dynamic structure and motion recovery as well as matching tensors for images of dynamic features. General, as well as highly constrained, dynamic scenarios, involving monocular or stereo views, have been investigated.

In this paper, we consider that the observed scene has both a static and a dynamic part. The static part is unconstrained (but has to be 3D) whereas on the other hand, as in [START_REF] Avidan | Trajectory triangulation: 3D reconstruction of moving points from a monocular image sequence[END_REF][START_REF] Han | Reconstruction of a scene with multiple linearly moving objects[END_REF][START_REF] Han | Multiple motion scene reconstruction from uncalibrated views[END_REF][START_REF] Sturm | Structure and motion for dynamic scenes -the case of points moving in planes[END_REF][START_REF] Wolf | On projection matrices P k → P 2 , k = 3, • • • , 6, and their applications in computer vision[END_REF][START_REF] Wolf | Join tensors: On 3D-to-3D alignment of dynamic sets[END_REF], we consider that dynamic features move along straight lines, termed motion lines. To further constrain the scenario, we consider that all motion lines lie on a motion plane and converge to an incidence point. Figure 2 illustrates this setting. Note that no assumption is made about the camera motion, which rules out background subtraction techniques, and that the camera is not assumed to be calibrated. A real-world instance of this scenario may be the motion of points arising from roadways seen from above, as for instance, by a moving surveillance video camera, see figure 1.

Figure 1: Selected frames of a video sequence consisting of 60 frames. Modeling the geometry of images provided by a moving camera observing a roadway is one application of our scenario. The difficulties of handling such sequences come from the fact that the scene consists of both static and dynamic features and that there is no constraint on the unknown camera motion.

This scenario fits into less constrained cases previously examined [START_REF] Avidan | Trajectory triangulation: 3D reconstruction of moving points from a monocular image sequence[END_REF][START_REF] Han | Reconstruction of a scene with multiple linearly moving objects[END_REF][START_REF] Han | Multiple motion scene reconstruction from uncalibrated views[END_REF][START_REF] Sturm | Structure and motion for dynamic scenes -the case of points moving in planes[END_REF][START_REF] Wolf | On projection matrices P k → P 2 , k = 3, • • • , 6, and their applications in computer vision[END_REF][START_REF] Wolf | Join tensors: On 3D-to-3D alignment of dynamic sets[END_REF]. The corresponding dynamic structure-and-motion algorithms and matching tensors may therefore be used. The main drawback is that they require, in general, a number of point correspondences that may not be well-adapted for e.g. robust estimation based on random sampling techniques, that is most of the time required to devise practical systems.

Moreover, the following drawbacks arise. The method proposed in [START_REF] Avidan | Trajectory triangulation: 3D reconstruction of moving points from a monocular image sequence[END_REF] requires that the camera motion is known and that point correspondences over 5 images are provided. The solutions proposed in [START_REF] Han | Reconstruction of a scene with multiple linearly moving objects[END_REF][START_REF] Han | Multiple motion scene reconstruction from uncalibrated views[END_REF] rely on the fact that observed features have constant velocity, as well as most applications, apart from segmentation tensors, provided in [START_REF] Wolf | On projection matrices P k → P 2 , k = 3, • • • , 6, and their applications in computer vision[END_REF]. In [START_REF] Sturm | Structure and motion for dynamic scenes -the case of points moving in planes[END_REF][START_REF] Wolf | Join tensors: On 3D-to-3D alignment of dynamic sets[END_REF], 3D views of the scene are required, which implies the use of two or more synchronized cameras. The H-tensor of [START_REF] Shashua | Homography tensors: On algebraic entities that represent three views of static or moving planar points[END_REF] necessitates at least 3 images of static or dynamic points to be computed.

We show that much simpler matching tensors and dynamic structure and motion algorithms may be derived for the case studied in this paper.

Firstly, in §3, we examine the purely dynamic case, i.e. when only dynamic points are observed and when camera motion is unknown. We show that in the two-view case, dynamic structure and motion may be described by a fundamental matrix-like tensor that we call the C-tensor. Standard techniques, such as robust estimation [START_REF] Torr | Motion Segmentation and Outlier Detection[END_REF] and maximum likelihood estimation through bundle adjustment [START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF] can then be applied in a straightforward manner to recover this tensor. We then show how dynamic motion in the multiple-view case can be modeled using a network of constrained C-tensors. A means to consistently estimate this geometry is provided. Dynamic structure in this case is also examined.

Secondly, in §4, we investigate the links between the previously-derived dynamic structure and motion and the projective registration (i.e. static motion) of the images. We give means for constrained estimation of camera motion.

Experimental results on real images may be seen throughout the paper and in §5. The dynamic scenario considered in this paper -unconstrained static points and points moving independently in a pencil of lines converging to the incidence point and lying on the motion plane. The motion plane and the incidence point are assumed to be static. However, most results derived in this paper are still valid without this assumption.

In most practical cases, the incidence point lies at infinity since motion lines are parallel.

Background and Notation

We consider sets of 3D points, each of them denoted as U , that may be split into dynamic and static points respectively denoted as X and Q. Corresponding time-varying 3D coordinates are respectively denoted by U, U , U , . . . , X, X , X , . . . and Q. Images of these points are respectively denoted by u, u , u , . . . , x, x , x , . . . and q, q , q , . . . . Figure 2 illustrates some of these notations. The incidence point is denoted by B. It has coordinates B and projects to b, b , b , . . . . It lies on the motion plane π that has coordinates π T ∼ ( π 1). The projective space of dimension d is denoted by P d . Everything is homogeneous (i.e. defined up to scale).

Purely Dynamic Views

Here we restrict to the case where only dynamic points can be observed from the scene. We assume that the different views are not registered, i.e. projection matrices are not available. We derive dynamic matching tensors for the two-then the multiple-view case.

Figure 3 shows a toy example overlaid with dynamic features.

Two Views: The 7-dof C-Tensor

Derivation. We propose a way to derive the C-tensor, encapsulating the dynamic twoview motion. Alternatively, other means could be used, such as P 3 → P 2 projection matrices within the framework of [START_REF] Wolf | On projection matrices P k → P 2 , k = 3, • • • , 6, and their applications in computer vision[END_REF] or similarly to the join tensors of [START_REF] Wolf | Join tensors: On 3D-to-3D alignment of dynamic sets[END_REF].

Let H be the unknown homography induced by the motion plane between the two views considered. Using H, we may predict the projection x of X in the second view if X was static as x ∼ Hx. The image line m of the motion line associated to X can then be obtained in the second view as the line joining b and x : m ∼ [b ] × Hx, where [b ] × is the (3 × 3) skew-symmetric cross-product matrix. Obviously, a necessary condition for X and X to be instances of the same dynamic point X is that x lies on m , which yields:

x T Cx with C ∼ [b ] × H, (1)
where we call C the 7-dof C-tensor (see below). It encapsulates the image signature of the dynamic two-view motion for the scenario previously described.

It is straightforward to see that the C-tensor has the same algebraic structure as the fundamental matrix. More precisely, the following analogy may be established. The projections of the incidence point play the roles of the epipoles while the 1D homography between the two motion line pencils corresponds to the epipolar transformation.

Properties. From the above-proposed analogy, several properties of the C-tensor may be easily derived. The C-tensor is rank-2 and has 7 dof. The projection of the indidence point in the first image, respectively in the second image, is the right null-space, respectively the left null-space, of the C-tensor: Cb = C T b = 0 (3×1) . The extended motion line transformation G (a 5-dof 2D line-to-line homography relating the motion line pencils) is linked to the C-tensor as: Similarly, a canonic plane homography, denoted as H can be recovered as well as a 3-dimensional set of 2D homographies H a consistent with C:

G ∼ C [b] × and C ∼ G [b] × . (2)
H ∼ [b ] × C and H a ∼ H + b a T . (3)
Note that H a is a set of 2D homographies containing the plane homography H induced by the motion plane. We will see in §4.1 that when the fundamental matrix F (weak calibration of the cameras) is available, it is possible to recover the unknown H by computing the intersection of the family H a and the 3-dimensional family of plane homographies defined by F.

Estimation. Another consequence of the analogy between the C-tensor and the fundamental matrix is that one can apply any two-view projective structure and motion algorithm to estimate C. For instance, we use the 8 point algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] embedded in a RANSAC-based robust estimation scheme [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis[END_REF] to compute an initial guess of C, that we further refine using uncalibrated two-view bundle adjustment. In this case, the projective depths of points represent in fact their displacement along the motion lines. Figure 4 shows the result of computing pair-wise C-tensors.

Figure 4: Motion line pencils estimated using two-view projective bundle adjustment. For the middle image, we compute the motion line pencils with respect to two C-tensors (with the first and with the third images). Note the significant discrepancy between them. This discrepancy will be eliminated by a consistent parameterization of the multiple-view case, see §3.2. Note also that points lying on the overtaking car have been discarded as outiers.

Multiple Views: The 5-dof C-Tensor

We show that the relationships between n unregistered dynamic views are contained in n -1 C-tensors submitted to some additional consistency constraints, which express the fact that the incidence point is unique, regardless of the current time instant.

Degrees-of-freedom analysis and derivation. From the previous section, we already know that for n = 2 views, the dynamic geometry has 7 dof and is represented by a 7dof C-tensor, say C. Consider now a third view of the same scene which shares dynamic features with at least one of the two other views, say the second one. One can compute the C-tensor C between the second and the third view. However, one has to remember that for a given time instant, the incidence point has a fixed position. Therefore, C and C have to share the second image b of the incidence point for being consistent, which provides 2 constraints and leaves 7 + 7 -2 = 12 dof for the dynamic geometry for n = 3.

It is then straightforward to derive that the n-view case has 7 + 5(n -2) = 5n -3 dof.

The dynamic geometry of a set of n images can therefore be conveniently modeled by a 7-dof C-tensor between two reference views and a network of n -2 5-dof C-tensors. A 5-dof C-tensor is a C-tensor with one constrained incidence point. A means to compute such a constrained C-tensor, given its right kernel, is provided in the next paragraph.

Concerning the dynamic structure, each point has 2 + (n -1) = n + 1 dof corresponding to its position in the motion plane and n -1 motions along its motion line.

Threading C-tensors. Once a solution has been obtained for the 7-dof C-tensor modeling the dynamic geometry of two particular images, subsequent 5-dof C-tensors have to be computed given one projection of the incidence point, e.g. with their left or right kernel known. Enforcing these consistency constraints when threading C-tensors is important since only 5 point correspondences instead of 7 are necessary to solve for the constrained tensor, as shown below. Moreover, the solution obtained is consistent and may be refined directly without prior correction using non-linear methods to obtain, e.g. a maximum likelihood solution.

We propose a linear algorithm inspired by [START_REF] Hartley | Minimizing algebraic error[END_REF] for estimating 5-dof C-tensors using 5 or more point correspondences while enforcing the consistency constraints. Using equation ( 1) and the factorization (2) of the C-tensor, we may write: ). An efficient solution to obtain this decomposition is given in [START_REF] Hartley | Minimizing algebraic error[END_REF]. By replacing into equation ( 4), we obtain: x T Ḡȳ = 0 with ( Ḡ(3×2) 0 (3×1) ) = GUΣ and (ȳ T (2×1) y) = V T x. Note that Ḡ is defined by 6 homogeneous parameters ḡ ∈ P 5 , which is consistent with the fact that G has 5 dof1 . Expanding this equation leads to an homogeneous linear system for ḡ. Note that 5 equations (hence 5 point correspondences) are sufficient to solve for ḡ in the least squares sense. From Ḡ, one can further obtain C as: C ∼ ( Ḡ(3×2) 0 (3×1) )V T .

x T G [b] × x = 0, (4)
Figure 5: Motion line pencils estimated using two-view projective bundle adjustment between the first and the second views and constrained estimation of a 5-dof C-tensor between the second and third views. Note that the motion lines pencils are perfectly aligned in the second image.

Estimation. The 7-dof C-tensor between two particular views may be estimated using the 8 point algorithm. Other 5-dof C-tensors may be computed using the above-described constrained method, possibly embedded in a RANSAC-based robust estimation process [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis[END_REF]. Figure 5 shows the result of such a constrained computation. This provides an initial guess of the dynamic registration, which may be refined as follows.

As in the case of multiple-view bundle adjustment, minimizing the reprojection error, i.e. the discrepancy between measured and predicted features yields the maximum likelihood estimator. We employed such a means directly for the two-view case since there was a direct analogy. For the multiple-view case however, one can not directly use standard bundle adjustment techniques for the following reasons. Firstly, the projective multipleview motion has 11n -15 dof whereas dynamic motion involves 5n -3 dof. Secondly, a reconstructed static point has 3 dof whereas a dynamic point has n + 1 dof (provided it is visible in every image considered). Therefore, the problem must be specifically parameterized.

Mixing Static and Dynamic Features

In this section, we consider that enough static points may be used to perform a weak calibration of the cameras, i.e. to projectively register them. A classical means for such a registration is to compute static structure and motion between two particular views and iteratively register the other views using 3D-to-2D point correspondences. We derive the links between the dynamic structure and motion given in §3 and classical static structure and motion and show how the above-mentioned registration algorithm can be constrained by known C-tensors.

Two Views

We represent the projective two-view motion by the fundamental matrix F. Firstly, we investigate the link between F and the C-tensor C. Secondly, we show that if F and C are known, the plane homography H induced by the motion plane may be recovered and we give a closed-form solution in terms of F and C as well as a means to use standard homography estimation algorithms.

The Link Between C and F

To establish this link, we consider the plane homography H induced by the motion plane. This homography can be written in terms of the C-tensor as H a , see equation [START_REF] Costeira | A multi-body factorization method for motion analysis[END_REF], where the unknown 3-vector a is used to span the space of 2D homographies consistent with C. The fundamental matrix can be formed from any plane homography as F ∼ [e ] × H and in particular H a which yields:

F ∼ [e ] × ([b ] × C + b a T ). (5)
This equation, that we call the F-C-consistency constraint, shows that F has only 5 dof corresponding to the right epipole and the equation a of the motion plane. It is equivalent to the fact that F and C share a 2D homography. Therefore, given C, 5 point correspondences should be enough to estimate F. However, due to the non-linearity of equation ( 5) for the unknowns e and a, we can not estimate F linearly using 5 point correspondences. Another solution is to use the fact that the incidence point B is a static point and that therefore, b and b give one constraint on F through the fundamental equation b T Fb = 0. A minimum of 6 other static points are then sufficient to estimate F.

Retrieving H

Given the C-tensor and the epipolar geometry, it is possible to recover the plane homography induced by the motion plane π. The following two paragraphs give respectively a closed-form solution in terms of the C-tensor and the fundamental matrix and a more physically meaningful solution taking feature positions into account.

A closed-form solution. Let x, x be the projections of a 3D point X ∈ π. Obviously, x lies at the intersection of the motion line and the epipolar line associated to x:

x ∼ (Cx) × (Fx).

We claim that H can be recovered as:

H ∼ (C × F) • diag (C × F) -1 i c i × i f i , (7)
where C × F is the column-wise cross-product of C and F. The proof of this result is omitted and may be found in [START_REF] Bartoli | The geometry of dynamic scenes -on coplanar and convergen linear motions embedded in 3D static scenes[END_REF].

Figure 6: Dynamic points transferred between images using plane homographies associated with the motion plane, i.e. as if they were static. These plane homographies have been estimated by hallucinating point correspondences (see text) and minimizing the symmetric squared distance between measured and transferred features. The final error is 0.35 pixels while the initial guess provided by the closed-form solution (not shown here) has an error of 2.77 pixels.

A physically meaningful solution. More physically meaningful means to estimate H, but more computationally expensive, can be obtained by hallucinating static point correspondences that lie on the plane using equation [START_REF] Han | Reconstruction of a scene with multiple linearly moving objects[END_REF]. Any standard method can then be used to solve for H by minimizing a given criterion, see e.g. [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]:

1. hallucinate point correspondences as:

(x, (Cx) × (Fx)) , x , (C T x ) × (F T x )
2. use any standard method to estimate H.

For example, we have chosen to non-linearly optimize the following cost function using the Levenberg-Marquardt algorithm initialized by the previously-given closed-form solution:

min H x↔x d 2 (Hx, (Cx) × (Fx)) + d 2 H -1 x , (C T x ) × (F T x ) .
Experimental results can be seen on figure 6.

Multiple Views

In the registered case, 2 views entirely fix the dynamic motion, since the incidence point and the motion plane may be determined uniquely. On the other hand, each additional view of a dynamic point adds 1 dof, corresponding to its position on its motion line at the time instant the picture was taken, as in the purely dynamic case. More details are given in [START_REF] Bartoli | The geometry of dynamic scenes -on coplanar and convergen linear motions embedded in 3D static scenes[END_REF], in particular, the link between C and P via the P-C-consistency constraint and the consistent estimation of these entities are investigated.

Experimental Results Using Real Images

We compute dynamic and static structure and motion on a 60-frame sequence from which sample images are shown on figure 1. We select dynamic and static features by hand on the first image and automatically track them through the sequence using a correlationbased technique. We then used key frames 0, 10, 20, . . . shown on figure 1.

We first perform dynamic structure and motion by sequentially computing C-tensors as described in §3.2. We then perform constrained static structure and motion as described in §4. Lastly, we use these results to recover the plane homographies associated to the motion plane between key frames, as indicated in §4.1.2. Such homographies allow to transfer dynamic features and predict their position in another camera position and another time instant as if they were static. The result of such transfers is shown on figure 7.

Figure 7: Plane homographies recovered using dynamic and static structure and motion allow to transfer vehicle positions from one key frame to the others as if they were static. The last two images are zooms on frames 0 and 10 respectively. The first one shows manually clicked static points (lying on the motion plane) while the second one shows the transfer of these points using the recovered plane homography (computed using dynamic points only).

Figure 7 also shows that the first homography (i.e. between frames 0 and 10) is relatively accurate since static point positions after transfer seem visually good.

The main problem that we encountered was the computation of the initial 7-dof Ctensor between the two first frames. Indeed, one may observe that all vehicles have roughly the same speed, which therefore induces a point-to-point homography between dynamic features of these two frames. There was therefore a 2-dof ambiguity on the computation of the 7-dof C-tensor. Instead, we computed a 5-dof C-tensor constrained by the projection of the incidence point in the first image. This projection was obtained by intersecting support lines of white bands on the ground.

Conclusion

We addressed the case of a specific dynamic scenario describing the motion of point features along lines converging to the same point and lying onto a motion plane. We show that very simple matching tensors that we call C-tensors, similar to fundamental matrices, exist. We show how to constrain static structure and motion by its dynamic counterpart. Plane homographies associated with the motion plane can then be recovered from dynamic features only. Experimental results show that this approach is feasible in practice and may be used to model e.g. surveillance video cameras observing roadways.

We believe that these geometrical features may be successfully used to devise completely automatic systems for vehicule tracking and camera motion estimation. Among issues for further work, self-calibration of the camera by considering that in practice the incidence lies most of the time at infinity, could be examined.

More details may be found in [START_REF] Bartoli | The geometry of dynamic scenes -on coplanar and convergen linear motions embedded in 3D static scenes[END_REF].

  Figure2: The dynamic scenario considered in this paper -unconstrained static points and points moving independently in a pencil of lines converging to the incidence point and lying on the motion plane. The motion plane and the incidence point are assumed to be static. However, most results derived in this paper are still valid without this assumption. In most practical cases, the incidence point lies at infinity since motion lines are parallel.

Figure 3 :

 3 Figure 3: Dynamic points used for the experiments on the toy images. Note that four points are lying on a car which overtakes another one in the second image, and therefore do not fulfill the dynamic motion associated to the other points.

  To understand the above expressions, consider a motion line m in the first image: [b] × m is a point on this line (b is interpreted as a line that does not contain the point b) and C [b] × is the corresponding motion line in the other image. A similar reasoning may be done to understand the expression of C from G.

  where G contains the unknown motion line pencil transformation and b is the known projection of the incidence point in the first image. Let us see how to solve for G. Let [b] × ∼ UΣV T be a singular value decomposition of [b] × where Σ = diag(1, σ, 0

ḡ is the row-wise vectorization of Ḡ.