Bernard Bonnard
email: bernard.bonnard@u-bourgogne.fr

Alain Jacquemard
email: dbacalain.jacquemard@u-bourgogne.fr

Mohab Safey
email: mohab.safey@lip6.fr

El Din

Thibaut Verron
email: thibaut.verron@lip6.fr

Determinantal Sets, Singularities and Application to Optimal Control in Medical Imagery

Keywords: Polynomial system solving, Real geometry, Applications

Control theory has recently been involved in the field of nuclear magnetic resonance imagery. The goal is to control the magnetic field optimally in order to improve the contrast between two biological matters on the pictures. Geometric optimal control leads us here to analyze meromorphic vector fields depending upon physical parameters, and having their singularities defined by a determinantal variety. The involved matrix has polynomial entries with respect to both the state variables and the parameters. Taking into account the physical constraints of the problem, one needs to classify, with respect to the parameters, the number of real singularities lying in some prescribed semialgebraic set. We develop a dedicated algorithm for real root classification of the singularities of the rank defects of a polynomial matrix, cut with a given semi-algebraic set. The algorithm works under some genericity assumptions which are easy to check. These assumptions are not so restrictive and are satisfied in the aforementioned application. As more general strategies for real root classification do, our algorithm needs to compute the critical loci of some maps, intersections with the boundary of the semi-algebraic domain, etc. In order to compute these objects, the determinantal structure is exploited through a stratification by the rank of the polynomial matrix. This speeds up the computations by a factor 100. Furthermore, our implementation is able to solve the application in medical imagery, which was out of reach of more general algorithms for real root classification. For instance, computational results show that the contrast problem where one of the matters is water is partitioned into three distinct classes.

INTRODUCTION

Motivations and problem description. Nuclear Magnetic Resonance (NMR) is a powerful tool in medical imagery. In order to distinguish two biological matters on a picture, it is required to optimize the contrast between the two matters. Because of its importance in medical sciences, this contrast imaging problem has received a lot of attention. The pioneering work of [3] has established geometric optimal control techniques as a major tool for designing optimal control strategies for the problem of improving the contrast.

These strategies depend on the biological matters under study. In NMR imagery the main physical parameters involved are the longitudinal and transversal relaxation times of each matter. This approach is formalized in [2]. It requires us to solve the following real root classification problem.

Consider a k × k matrix M whose coefficients are polynomials in Q[X 1 , . . . , X n , G 1 , . . . , G t ], and assume that n = (kr + 1) 2 with r ∈ {1, . . . , k -1}. Let π : C n ×C t → C t be the canonical projection. Let V r ⊂ C n × C t be the set of points at which M has rank r, and let V be the union of the singular locus of V r and of the critical points of π restricted to V r . Generically, this variety has dimension n+t -(k -r +1) 2 = t. Also consider a semi-algebraic set B in R n × R t . Assume that B has non-empty interior and that there exists a Zariski-open set O ⊂ C t such that V ∩ π -1 (g) is a non-empty finite set for g ∈ O. Further assume that the projection π restricted to V ∩ B is proper ([10, Def. 2.10.1]). We aim at describing connected open sets C 1 , . . . ,C ⊂ R t such that i=1 C i is dense in R t (for the Euclidean topology) and the cardinality of V ∩ π -1 (g) inside B ∩ π -1 (g) is invariant when g ranges over C i .

For our application, the size of the matrix is 4, and the target rank is k -1 = 3. The number of variables and parameters are n = 4 and t = 4 for the general case. Since the system is homogeneous in the parameters, we may set one of the parameters to 1, reducing the problem to t = 3. An important particular case is when one of the matters is water; then the corresponding relaxation times are 1, and the number of parameters is t = 2.

Such a determinantal structure is general enough to design dedicated algorithms. We also mention that the optimal control problems in [6] lead to algebraic classification problems with similar structures.

State-of-the-art. The modeling through an optimal control problem is introduced in [2]. The so-called Bloch modeling and saturation method for tackling this problem is developed therein.

In [2], four experimental important cases are studied (all parameters of the classification problem are fixed). Among other properties, it has been observed there that the number of singularities is constant when water is involved. This led to the following questions: 1. Is this number of singularities preserved for any choice of a second matter, the first one being water?

2. If not, how many different classes of pairs of matters can we distinguish through the analysis of those singularities? Answering these questions leads to the real root classification problem described above. Symbolic computation techniques are good candidates to solve them.

Properties of Cylindrical Algebraic Decomposition (CAD) adapted to a given polynomial family allow to solve real root classification problems. Hence the CAD algorithm [START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition[END_REF] can be used in our context. However, the complexity of computing a CAD is doubly exponential in the number of variables ( [4,9]); its implementations are usually limited to non-trivial problems involving 4 variables and cannot tackle our application.

The complexity of computing a CAD can be much improved when taking into account equational constraints (see e.g. [START_REF] Mccallum | On projection in CAD-based quantifier elimination with equational constraint[END_REF]). In the context of real root classification problems, this leads us to take advantage of the presence of equations to compute closed sets in the parameter space (R t using our notation) containing the boundaries of the regions C 1 , . . . ,C , hence substituting the recursive (doubly exponential) projection steps of CAD with more involved projection techniques. In the past ten years, several works have focused on this problem [START_REF] Lazard | Solving parametric polynomial systems[END_REF][START_REF] Yang | A complete algorithm for automated discovering of a class of inequality-type theorems[END_REF] using various computer algebra tools such as Gröbner bases, regular chains, etc. We also mention [START_REF] Rodriguez | Data-discriminants of likelihood equations[END_REF] which uses evaluation/interpolation techniques to compute those closed sets in the parameters space.

While the implementation of [START_REF] Yang | A complete algorithm for automated discovering of a class of inequality-type theorems[END_REF] is able to solve our classification problem for the case of water, none of the implementations were able to classify the singular locus of V in the general case (the number of parameters is 3).

Our strategy is not as general as the ones in [START_REF] Yang | A complete algorithm for automated discovering of a class of inequality-type theorems[END_REF] or [START_REF] Lazard | Solving parametric polynomial systems[END_REF]. It exploits properties of sets defined by minors of matrices with polynomial entries. Such structures have been used for computing sample points in each connected component of the real trace of determinantal varieties [START_REF] Henrion | Real root finding for low rank linear matrices[END_REF][START_REF] Henrion | Real root finding for rank defects in linear Hankel matrices[END_REF][START_REF] Henrion | Real root finding for determinants of linear matrices[END_REF] or for solving linear matrix inequalities [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF]. These works are based on dedicated strategies for computing critical loci of some projections restricted to determinantal varieties. Such computations are naturally related to real root classification problems and real quantifier elimination (see e.g. [START_REF] Hong | Variant quantifier elimination[END_REF]). Finally, our computations rely on Gröbner bases. Several works have shown some connection between Gröbner bases and determinantal ideals [START_REF] Faugère | On the complexity of the Generalized MinRank Problem[END_REF] and critical point computation [START_REF] Faugère | Critical points and Gröbner bases: the unmixed case[END_REF][START_REF] Spaenlehauer | On the complexity of computing critical points with Gröbner bases[END_REF]. Main results. Our main results are twofold:

• an algorithm solving the special real root classification problem described above and which exploits the determinantal structure of the input data arising in contrast imaging problems;

• its successful use for solving the challenging application to the contrast problem in the general case;

• answers to the questions raised by the experimental data involving water: the answer to question 1. is no, and the answer to question 2. is that there are 3 classes of second matters that we can separate, depending on whether there are 1, 2 or 3 singularities.

We start by describing our algorithmic contribution. Recall that we are given a matrix, denoted by M, with polynomial entries. As in [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF][START_REF] Henrion | Real root finding for low rank linear matrices[END_REF][START_REF] Henrion | Real root finding for rank defects in linear Hankel matrices[END_REF][START_REF] Henrion | Real root finding for determinants of linear matrices[END_REF], it is based on splitting computations according to the rank of M. More precisely, in order to solve our real root classification problem, we need to identify where the number of real solutions inside B of the determinantal system describing V changes depending on the values of the parameters. The real root classification problem we want to solve involves inequalities defining a semi-algebraic set with non-empty interior. In this context, we use standard tools from real geometry, such as Thom's first isotopy lemma, which reduce our classification problem to computing the singular points of V , the critical points of the projection of the parameter space restricted to V , and the intersection of V with the boundary of the semi-algebraic set B.

This computation may be difficult because generically, the variety V r has singularities corresponding to points where rank(M) < r: this is proved using Bertini's theorem and [5, Prop. 1.1], as in [START_REF] Henrion | Real root finding for low rank linear matrices[END_REF]Prop. 2]). Hence, observe that the variety V is naturally split according to the rank of M. This is the very basic idea on which our algorithm relies: we compute critical loci of the projection on the parameter space restricted to the variety V by distinguishing those points at which M has rank less than r from those at which M has rank exactly r. Incidentally, it raises the question of how these higher rank deficiencies should be interpreted from the application viewpoint. To the best of our knowledge, this question was never raised in the optimal control community.

Our algorithms need to compute projection of algebraic sets, which is done using elimination algorithms such as Gröbner bases or triangular sets for example. We have performed experiments for both these tools, using the package FGb [START_REF] Faugère | FGb: A Library for Computing Gröbner Bases[END_REF] in Maple and an implementation of F5 [START_REF] Faugère | A new efficient algorithm for computing Gröbner bases without reduction to zero (F 5 )[END_REF] for Gröbner bases, and using the package RegularChains [START_REF] Lemaire | The regularchains library[END_REF] in Maple for triangular sets.

Regarding the contrast imaging problem, we illustrate the behaviour of our algorithm in the case of water, giving the whole classification. Using Gröbner bases to perform the eliminations, the computation takes 10 s on an 2 GHz Intel Xeon CPU. The Re-alRootClassification command of the Maple RegularChain library needs 1600 s to find this classification.

We also ran our algorithm on the general case. While none of the available implementation is able to tackle this classification problem directly, ours can find the polynomials separating the open sets C i within 4 h using FGb, or 2 min using F5, and the projection step of the CAD can be done in 4 h. We see similar speed-ups when using triangular sets to perform the elimination.

This illustrates how our dedicated algorithms take advantage of the special structure of the problem, to achieve speed-ups when compared with more general techniques.

Conclusion.

We propose an algorithmic strategy refining general roots classification strategies for the case of singularities of a determinantal variety, under genericity hypotheses. We give an overview of the results for the application at the end of the paper1 . We were able to give a full classification in the case of water, answering the questions raised by the experiments. For the general case, the separating polynomials were found, it remains to perform the analysis of the subdivisions in order to obtain the full classification. This work also raised questions concerning the interpretation of higher rank deficiencies from the viewpoint of control theory.

Structure of the paper. In Section 2, we present the mathematical background around NMR imagery and the contrast problem. Sec-tion 3 deals with the dedicated classification algorithm. Finally, in Section 4, we report on experimental results obtained when solving the application to the contrast imaging problem.

We acknowledge the support from the ANR-DFG research program Explosys (Grant No.ANR-14-CE35-0013-01; GL203/9-1). Mohab Safey El Din is member of and supported by Institut Universitaire de France.

MODELING THE DYNAMICS

The model we describe below has been introduced in [3] in order to apply techniques from geometric optimal control theory to the control of the spin dynamics by NMR. Up to some normalization, each spin 1/2 particle is governed by the Bloch equation

     ẋ = -Γx + u y z ẏ = -Γy -u x z ż = γ(1 -z) + u x y -u y x,
where the state variable q = (x, y, z) represents the magnetization vector which must lie in the Bloch ball defined by |q| ≤ 1, and the parameters (Γ, γ) are related to the physical relaxation times. The parameters must also satisfy 2Γ ≥ γ > 0. The control u = (u x , u y ) represents the magnetic field whose magnitude is bounded by a maximum value µ.

In the context of the contrast imaging problem, this leads to the simultaneous control of two non-interacting spins with different relaxation time parameters. The contrast by saturation method consists in bringing the magnetization vector of the first spin toward the center of the Bloch ball while maximizing the modulus of the magnetization vector of the other matter. The matter with a zero magnetization is black on the picture, while the other matter with a maximum modulus of the magnetization vector is bright.

Using the symmetry of revolution [2] which allows to eliminate one state variable for each matter, we obtain the system

     ẏ1 = -Γ 1 y 1 -u x z 1 ż1 = γ 1 (1 -z 1 ) + u x y 1 ẏ2 = -Γ 2 y 2 -u x z 2 ż2 = γ 2 (1 -z 2 ) + u x y 2 , ≤ µ (1)
and the optimal control problem is: starting from the equilibrium point N = ((0, 1), (0, 1)), saturate the first spin, that is q 1 (T ) = 0, where T is the transfer time while maximizing |q 2 (T )| 2 , where |q 2 (T )| represents the final contrast. It is a standard Mayer problem in optimal control, studied in [3] through the analysis of the Hamiltonian dynamics given by the Pontryagin Maximum Principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF].

We summarize this analysis below. Writing (1) as q = F(q) + u G(q), |u| ≤ µ, the optimality conditions associated with the Maximum Principle lead us to construct the optimal solution as a concatenation of bang-arcs where the control is u = ±µ, and singular arcs solutions of X e = F + u s G where the control u s is the rational fraction -D /D with

D = det(F, G, [G, F], [[G, F], G]) D = det(F, G, [G, F], [[G, F], F]),
where [ , ] denotes the Lie bracket of vector fields. Explicitly, with

d i = γ i -Γ i (i ∈ {1, 2}), D = det -Γ 1 y 1 -z 1 -1 d 1 z 1 -Γ 1 2 d 1 y 1 -γ 1 z 1 y 1 d 1 y 1 -2 d 1 z 1 +Γ 1 -d 1 -Γ 2 y 2 -z 2 -1 d 2 z 2 -Γ 2 2 d 2 y 2 -γ 2 z 2 y 2 d 2 y 2 -2 d 2 z 2 +Γ 2 -d 2 .
The localization of the singularities of {D = 0} inside the Bloch ball is important to understand the geometry of the hypersurface, as well as the dynamics of the vector field X e which is closely linked to the presence of such singularities. Indeed, generically when approaching the surface D = 0 along a singular arc, the con-trol u = -D /D is such that |u| > µ. Hence the control policy switches from singular to bang. Therefore this surface is related to the complexity of the optimal law, as a concatenation of singular and bang arcs. See [3] for numerical simulations related to this phenomenon.

ALGORITHM

Classification strategy

We consider the polynomial algebra

Q[X, G] with variables X = (X 1 , . . . , X n ) and parameters G = (G 1 , . . . , G t ). Let F and H be families of polynomials in Q[X, G]. Let V R = V R (F), V = V C (F)
be respectively the set of zeroes of F in R n+t and in C n+t . Let B be the closed semi-algebraic set defined by H:

B = {(x, g) ∈ R n+t | ∀h ∈ H, h(x, g) ≤ 0},
and let B 0 = h∈H V C (h). Let π : C n+t → C t be the projection onto the affine space with coordinates G. Let sing(V ) be the singular locus of V , crit(π,V ) be the set of critical points of π restricted to V , and

K(π,V ) = π(sing(V ) ∪ crit(π,V )) ∩ R t .
Given a subset A of a real or complex affine space, A and ∂ A are used to denote the closure and the boundary of A for the Euclidean topology respectively.

Assume that the following hypotheses are satisfied:

H1 There exists a nonempty Zariski-open subset O 1 of C t such that for all g ∈ O 1 , the fiber V ∩ π -1 (g) is a nonempty finite subset of C n+t ;
H2 The restriction of the projection π to B is proper ([10, Def. 2.10.1]);

H3 The intersection V ∩ B 0 has dimension at most t -1 in C n+t ;
H4 The variety V is equidimensional of dimension t.

We want to find a non-zero polynomial

P ∈ Q[G] such that, on each connected component U of R t \ V R (P), for g ∈ U, the cardi- nality of V ∩ B ∩ π -1 (g) does not depend on g.
In Lemma 1, we describe a well-known strategy for computing these objects (see for example [START_REF] Lazard | Solving parametric polynomial systems[END_REF][START_REF] Yang | A complete algorithm for automated discovering of a class of inequality-type theorems[END_REF]). LEMMA 1. Let F and H be polynomial systems satisfying hypotheses H1, H2, H3 and H4. Let C B = π(V ∩B 0 ), U a non-empty connected open subset of R t which does not meet C B ∪K(π,V ), and g ∈ U. Then V ∩ π -1 (g) is finite, and

∀ g ∈ U, # V ∩ π -1 (g ) = # V ∩ π -1 (g) .
PROOF. We will construct a Whitney stratification of V ∩ B ([1, Def. 9.7.1]) with certain properties. First note that since V is tequidimensional by H4, V has real dimension at most t ([1, Prop. 2.8.2]). Let S =t be the intersection of the points where V R has local dimension t and of the interior of B. There exists a Whitney stratification (S i ) of the semi-algebraic set V ∩ B such that S =t is the union of strata of dimension t ([1, Th. 9.7.11]). Let S <t be the union of the other strata, they all have real dimension less than t. By construction, this is a semi-algebraic set which is the union of (V ∩ ∂ B) ⊂ (V ∩ B 0 ) and of the singular locus of V ∩ B, and it has dimension less than t (by H3 for V ∩ B 0 ). Its image through π has dimension less than t, and so it has codimension at least 1. Now consider S =t . By Hyp. H1, for any

g ∈ O 1 , π -1 (g) ∩ V is non-empty, hence π(V ) contains the non-empty Zariski-open set O 1 ⊂ C t . The intersection O 1 ∩ R t is a non-empty Zariski-open set of R t , contained in π(V ), hence π(V ) ∩ R t has real dimension t.
Let U 0 be its interior. The subset S =t ∩ π -1 (U 0 ) is a locally closed semi-algebraic set. If it is empty, then there is nothing to prove. Otherwise, by construction it has dimension t; and the projection π restricted to this subspace is proper, by Hyp. H2. Thom's isotopy lemma ( [START_REF] Coste | Thom's first isotopy lemma: a semialgebraic version, with uniform bound[END_REF]) states that for any nonempty connected open set U of R t not meeting K(π,V ), and for any g ∈ U, there exists a semi-algebraic diffeomorphism h = (h 0 , π) :

V ∩ B ∩ π -1 (U) ∼ -→ π -1 (g) ×U.
By Hyp. H1, if U is nonempty, π -1 (g) is finite, and the cardinality of the fibers is constant on U.

So computing the wanted decomposition of the parameter space is equivalent to computing a polynomial P ∈ Q[G] such that V (P) covers π(V ∩ B 0 ) and K(π,V ).

The determinantal problem

Let k be an integer greater than 1, r 0 ∈ {1, . . . , k -1}, and n = (kr 0 + 1) 2 . Let t ∈ N, and let M(X, G) be a k × k matrix with polynomial entries in n variables X = (X 1 , . . . , X n ) and t parameters G = (G 1 , . . . , G t ). As before, let π : C n+t → C t be the projection onto the affine space with coordinates G.

Let {h(X, G) ≤ 0 | h ∈ H} be a system of inequalities, with

H ⊂ Q[X, G]. Let V -1 = /
0 by convention, and for any r ∈ {0, . . . , k}, define the variety

V r = {(x, g) ∈ C n+t | rank(M(x, g)) ≤ r}
and the constructible set V =r = V r \V r-1 , that is the set of points at which the matrix M has rank exactly r.

Let V be the union of the singular locus of V r 0 and of the set of critical points of π restricted to V r 0 . We want to classify the cardinality of the real fibers by π of the semi-algebraic set

V ∩ {(x, g) | ∀h ∈ H, h(x, g) ≤ 0}.
Assume that V and H satisfy hypotheses H1, H2, H3 and H4. Further assume that:

H5 There exists a non-empty Zariski-open subset O 2 ⊂ C t such that V ∩ π -1 (O 2 ) = V r 0 -1 ∩ π -1 (O 2 );
H6 For any r ∈ {0, . . . , k -1}, the ideal defined by the (r + 1)minors of M is radical;

H7 For any r ∈ {0, . . . , k -1}, the variety V r is equidimensional with dimension n + t -(kr + 1) 2 . These properties are generic ([16, Prop. 2]). LEMMA 2. Assuming Hyp. H6, for any r ∈ {1, . . . , k}, V r-1 ⊂ sing(V r ).

PROOF. We will prove the following stronger statement: let g) ∈ V r-1 , then all partial derivatives of all (r + 1)-minors of M vanish at (x, g). This will prove that the Jacobian of V r at (x, g) is the zero matrix , and in particular has rank 0. By Hyp. H6 the ideal of all (r + 1)-minors of M is radical, so we can use the Jacobian criterion to characterize sing(V ), so (x, g) ∈ sing(V r ).

Let (x, g) ∈ V r-1 . For any (r + 1) × (r + 1)-submatrix of M, the result we want to prove depends only on the coefficients of the submatrix. So w.l.o.g., we may assume that r = k -1. Let (x, g) ∈ V k-2 , this means that all (k -1)-minors of M vanish at (x, g). Consider a matrix of polynomial indeterminates U:

M =   U 1,1 ... U 1,k . . . . . . U k,1 ... U k,k  
and let D be its determinant. Then for any i, j ∈ {1, . . . , k}, ∂ D ∂U i, j = (-1) i+ j • Mi, j (U) where Mi, j is the (k -1)-minor of M obtained by removing row i and column j.

For any i, j ∈ {1, . . . , k}, let m i, j (resp. M i, j ) be the coefficient at row i and column j of the matrix M (resp. the minor obtained by removing row i and column j from M). By the derivation chain rule, for any

v ∈ {X 1 , . . . , X n , G 1 , . . . , G t }, ∂ D ∂ v = ∑ k i, j=1 (-1) i+ j ∂ m i, j ∂ v • Mi, j (m), which equals ∑ k i, j=1 (-1) i+ j ∂ m i, j ∂ v • M i, j (X).
Since by hypothesis all (k -1) minors of M vanish at (x, g), all partial derivatives ∂ D ∂ v vanish at (x, g).

In the following subsections, we will describe two algorithms DeterminantCritVals and DeterminantBoundary, which, given such a matrix M, a target rank r 0 and inequalities H, compute a polynomial whose zeroes cover K(π,V ), and a polynomial whose zeroes cover π(V ∩ B 0 ) respectively. By Lemma 1, the zeroes of the product of these polynomials will subdivide the parameter space into connected components where the cardinality of real fibers is constant. These algorithms are probabilistic, because they will rely on the choice of generic linear forms to ensure linear independence. However, the algorithms could be made deterministic by testing that these linear forms are generic enough for our purpose, and repeating the random draw otherwise.

The algorithms will also need to compute the projection of algebraic sets onto coordinate subspaces. For this purpose, we assume that we are given a routine Elimination, which, given a system of polynomials

F ⊂ Q[V 1 , . . . ,V N ] and a set of variables V ⊂ {V 1 , . . . ,V N }, com- putes a system of generators of F ∩ Q[V ]
. Such a routine can be implemented using Gröbner bases or regular chains, for example.

Incidence varieties

We decompose the problem depending on the rank of the matrix. The classical technique that we use to model properties on the rank relies on incidence varieties. with the additional condition that the matrix (y i, j ) has rank kr.

The projection of V r onto the affine space with coordinates (X, G) is V r . Let (u 1,1 , . . . , u k-r,k ) ∈ C k(k-r) , we define the variety V r,u as the intersection of V r and the complex solutions of

u 1,1 ... u 1,k . . . . . . u k-r,1 ... u k-r,k • y 1,1 ... y 1,k-r . . . . . . y k,1 ... y k,k-r = Id k-r (3)
In the rest of Section 3, F r,u denotes the union of Equations ( 2) and (3). LEMMA 4. For any r ∈ {0, . . . , k -1}, the varieties V r and V r,u are birational ( [START_REF] Bochnak | Real algebraic geometry[END_REF]Sec. 3.4]).

PROOF. We need to define a morphism f : W → V r,u with W a non-empty Zariski-open subset of V r , and such that f is inverse to the projection V r,u → V r onto the affine space with coordinates (X, G). Let W be the open subset of V r defined as the non-vanishing locus of the top-left r-minor of M. Consider the block decompositions, where A, Y (1) and U (1) are r × r matrices:

M = A B C D Y = Y (1) Y (2) U = U (1) U (2)
. Over W , A is invertible, let ∆ = det(A). Let M/A be the Schur complement of A in M, Equations ( 2) and ( 3) can be rewritten

  ∆Id r A -1 B 0 M/A U (1) U (2)   • Y (1) Y (2) =   0 0 Id r  
We may restrict to the open subset of V r,u where Y (2) is invertible, then eliminating Y (1) yields that

Y (2) = (U (2) -U (1) A -1 B) -1 Y (1) = -1 ∆ A -1 BY (2)
which defines the wanted morphism W → V r,u .

PROPOSITION 5. Let r 1 = r 0 -1. Let ϕ : V r 1 ,u → C t be the projection onto the affine space with coordinates G. Assuming that hypotheses H1 to H7 hold, there exists a Zariski-open subset

U ⊂ C k(k-r 1 ) such that if u ∈ U ∩ Q k(k-r 1 ) , K(π,V r 1 ) = K(ϕ, V r 1 ,u ).
PROOF. Let P = (x, g, y) ∈ V r 1 ,u . If M(x, g) has rank less than r 1 , then by Lemma 2, (x, g) ∈ sing(V r 1 ), hence g ∈ K(π,V r 1 ).

Since M(x, g) has rank less than r 1 , its kernel L 1 has dimension at least kr 1 + 1. Equations ( 3) encode that the vectors y i given by the columns of matrix (y i, j ) generate a r 1 -dimensional vector space L 2 . So there exists y 0 ∈ L 1 ∩ L 2 , and for all a ∈ C, (x, g, y 1 + ay 0 , y 2 ) belongs to the fiber above (x, g) in V r 1 . So this fiber has dimension at least 1, while the generic fiber has dimension 0 by hypothesis H1. So (x, g) is a critical value of the projection of V r 1 onto R n+t , hence (g) ∈ K(ϕ, V r 1 ,u ).

So we may assume that M(x, g) has rank exactly r 1 . There is a r 1 × r 1 submatrix A of M(x, g) which is invertible, without loss of generality we may assume that it is the top-left r 1 ×r 1 submatrix. In an open neighborhood of (x, g), V =r 1 is described by the vanishing of the entries of M/A, that is the determinants of the (r 1 + 1) × (r 1 + 1) submatrices containing A. The same computations as in the proof of Lemma 4 give the following equations describing V r 1 ,u in the open neighborhood of (x, g, y) where ∆ = det(A) does not vanish:

     M/A = 0 Y (2) = (U (2) -U (1) A -1 B) -1 Y (1) = -1 ∆ A -1 BY (2) (4)
and the truncated Jacobian matrix in (X, Y) of this system can be written

  Jac X (M/A) 0 0 Id r 1 (k-r 1 ) 0 Id r 1 (k-r 1 )  
where Jac X (M/A) is the truncated Jacobian matrix in X of the

(k -r 1 ) 2 entries of M/A, which define V r 1 \{(x, g | ∆ = 0)} in C n+t .
By hypothesis H6, the ideal defined by the entries of M/A, which is a subideal of the ideal of all (r 1 + 1)-minors of M, is radical. Since the Schur complement appears by multiplication with invertible matrices with entries in the localized ring Q[X, g] ∆ (using the same notations as in the proof of Lemma 4):

Id r 1 0 -C Id k-r 1 • A -1 0 0 Id k-r 1 • M = Id r 1 A -1 B 0 M/A , Equations (4) describe the localization of F r 1 ,u in Q[X, g] ∆ ,
so this ideal is radical as well. So we can use the Jacobian criterion on V r 1 near (x, g) and on V r 1 ,u near (x, g, y). Both Jacobians matrices have the same rank and both varieties have the same local codimension (kr 1 ) 2 (by Lemma 4 and hypothesis H7), so

K(π,V =r 1 ) ∩ ϕ(V r 1 ,u ) = K(ϕ, V r 1 ,u ) ∩ π(V =r 1 )
The image ϕ(

V r 1 ,u ) ∩ π(V =r 1 ) is a Zariski-open subset O u of π(V =r 1 ). It remains to prove that if u is sufficiently generic, then all irreducible components of K(π,V =r 1 ) meet this open subset.
Let C 1 , . . . , C a be these irreducible components, and let (x 1 , g 1 ) ∈ π -1 (C 1 ), . . . , (x a , g a ) ∈ π -1 (C a ). For any (x, g) ∈ V =r 1 , the proof of [START_REF] Henrion | Real root finding for low rank linear matrices[END_REF]Prop. 4,Sec. 6] shows that there exists a non-empty Zariskiopen subset

U (x,g) ⊂ C k(k-r 1 ) such that if u ∈ U (x,g) ∩Q k(k-r 1 ) , then (x, g) ∈ O u ; namely, U is the set of u such that rank M(x,g) (u i, j ) = k.
Taking the finite intersection of the non-empty Zariski-open subsets U (x i ,g i ) for i ∈ {1, . . . , a} yields the wanted subset U.

Locus of rank exactly r 0

Recall that by H5, π(V ∩ V =r 0 ) has codimension at least 1, and that we want to compute a polynomial whose zeroes cover K(π,V ) and π(V ∩ B 0 ). So we may multiply the result by the equation of one hypersurface covering π(V ∩ V =r 0 ), it will naturally cover π(V ∩V =r 0 ) ∩ K(π,V ) and π(V

∩V =r 0 ) ∩ π(V ∩ B 0 ). Algorithm RankExactly: Input M ⊂ Q[X, G] k×k , r 0 ∈ {1, . . . , k -1} Output P 1 ∈ Q[G] \ {0} s.t. π(V ∩V =r 0 ) ⊂ V (P 1 ) Procedure PROPOSITION 6. Algorithm DeterminantCritVals is correct. PROOF. By definition, V ⊂ V r 0 . Using the decomposition V r 0 = V =r 0 ∪ V r 0 -1 , we decompose the variety V as V = (V ∩V =r 0 ) ∪ (V ∩V r 0 -1 ).
The subspace π(V ∩ V =r 0 ) is covered by the output of RankExactly, so we may restrict to V ∩ V r 0 -1 , which is the whole variety V r 0 -1 by Lemma 2. Recall that by hypothesis H7, V r 0 -1 is t-equidimensional.

By Prop. 5, in order to compute K(π,V r 0 -1 ), we can compute polynomials whose zeroes cover K(ϕ, V r 0 -1,u ) with u sufficiently generic instead.

By hypotheses H6, H7 and the proof of Prop. 5, V r 0 -1,u is tequidimensional and F r 0 -1,u is a set of generators of its ideal, so we can use the Jacobian criterion to compute equations defining K(ϕ, V r 0 -1,u ) .

Boundary

Algorithm DeterminantBoundary:

Input • M ⊂ Q[X, G] k×k • r 0 ∈ {1, . . . , k -1} • H ⊂ Q[X, G] (set of constraints on the variables) Output P b ∈ Q[G] \ {0} s.t. π(V ∩ B 0 ) ⊂ V (P b ) Procedure 1. res ← RankExactly(M, r 0 ) 2. Pick at random u 1 , . . . , u k(k-r 0 +1) = u ∈ Q k(k-r 0 +1) 3. F 0 ← F r 0 -1,u 4. For h in H do 5. F 1 ← F 0 ∪ {h} 6.
G ← Elimination(F, {X, Y}) 7.

Multiply res by 1 polynomial from G 8. End for 9. Return res PROPOSITION 7. Algorithm DeterminantBoundary is correct.

PROOF. As in Section 3.5, we write: V = (V ∩V =r 0 )∪(V ∩V r 0 -1 ). Since B 0 = h∈H V (h), the intersection V ∩ ∂ B is contained in the union of the varieties V ( F + h ) for h ranging over H, and the equation of the projections can be obtained with polynomial elimination. REMARK 8. For the real root classification problem, the subdivision is given by the product of the outputs of DeterminantCrit-Vals and DeterminantBoundary. In order to avoid repeating computations, we may skip the call to RankExactly in either subroutine (but not both), and initialize res to 1 instead.

THE CONTRAST PROBLEM

The case of water

With the notations of Section 2, the variety V is the complex algebraic variety defined by

D = ∂ D ∂ y 1 = ∂ D ∂ y 2 = ∂ D ∂ z 1 = ∂ D ∂ z 2 = 0.
With the notations of Section 3, we want to classify the singularities of the set of points where M has rank at most r 0 = 3. Our semialgebraic constraints are that the solutions are within the Bloch ball, that is

B : h 1 = y 2 1 + (z 1 + 1) 2 ≤ 1 h 2 = y 2 2 + (z 2 + 1) 2 ≤ 1.
Since the equations are homogeneous in Γ 1 , Γ 2 , γ 1 , γ 2 , and the parameters are supposed to be non-zero, we may normalize by setting γ 1 = 1. In the case where the first matter is water, we further simplify by setting Γ 1 = γ 1 = 1, leaving free the two parameters Γ 2 , γ 2 corresponding to the second matter. We recall that we also assume that 2 Γ 2 ≥ γ 2 and that (γ 2 , Γ 2 ) = (1, 1) = (γ 1 , Γ 1 ) (that is, the second matter is not water).

THEOREM 9. Consider the 9 polynomials:

f 1 = Γ 2 -1, f 2 = 3 Γ 2 -2 γ 2 -1, f 3 = 3 Γ 2 2 -5 Γ 2 γ 2 + γ 2 2 + 2 Γ 2 -2 γ 2 + 1, f 4 = 2 Γ 2 2 -5 Γ 2 γ 2 + 2 γ 2 2 -2 Γ 2 + 3 γ 2 , f 5 = 2γ 3 2 -(3Γ 2 + 11) γ 2 2 + 9Γ 2 + 6 -3Γ 2 2 γ 2 + 2 Γ 2 (Γ 2 + 2) (Γ 2 -1) , f 6 = Γ 2 -2 γ 2 + 1, f 7 = 2 Γ 2 -γ 2 -1, f 8 = γ 2 -2 + Γ 2 , f 9 = 2 Γ 2 2 -5 Γ 2 γ 2 + 2 γ 2 2 + 1.
The zeroes of their product divide the subset of R 2 defined by 2 Γ 2 > γ 2 > 0 into connected components where the cardinality of V R ∩ π -1 (γ 2 , Γ 2 ) is constant.

- 1 ∪ G + 1 (resp. G - 2 ∪ G + 2 )
. The circled numbers in each area correspond to the number of singularities in B for parameters in the area. Parameters in the red area are physically irrelevant.

Let ψ : (y 1 , z 1 , y 2 , z 2 ) → (-y 1 , z 1 , -y 2 , z 2 ) be the symmetry fixing Π = {y 1 = y 2 = 0}, and let us consider the semi-algebraic sets (see Fig. 1):

G - 1 = {0 < γ 2 < 2 Γ 2 , Γ 2 < 1, f 2 > 0}, G + 1 = {0 < γ 2 < 2 Γ 2 , Γ 2 > 1, f 2 < 0, f 4 > 0}, G - 2 = {0 < Γ 2 < 1, f 6 > 0, f 3 < 0}, G + 2 = {Γ 2 > 1, f 6 < 0, f 5 > 0}, G = G - 1 ∪ G + 1 ∪ G - 2 ∪ G + 2 .
THEOREM 10. For all (γ 2 , Γ 2 ) such that 2 Γ 2 > γ 2 > 0, the center O of the Bloch ball B is a singularity of {D = 0}. And provided (γ 2 , Γ 2 ) ∈ G, there exist at most two other singularities:

1. provided (γ 2 , Γ 2 ) ∈ G - 1 ∪ G + 1 there is one other singularity lying on Π ∩ B; 2. provided (γ 2 , Γ 2 ) ∈ G - 2 ∪ G + 2
there are two other singularities in B, ψ-symmetric, outside Π.

The configuration is illustrated in Figs. 1 and2. Observe that the number of singularities inside B is an invariant of the contrast problem. Two of the pairs of biological matters studied in [2], watercerebrospinal fluid (normalized parameters [γ 2 = 5 4 , Γ 2 = 25 3 ]) and water-fat (normalized parameters [γ 2 = 25 2 , Γ 2 = 25]) correspond to points outside G, and their invariant is 1 in both cases (see Fig. 2). But our results give answers to our guiding questions: there exist pairs of matters for which this algebraic invariant can differ; and any pair (water, matter) belongs to one of 3 classes, depending on whether the number of singularities inside B is 1, 2 or 3.

PROOF OF THEOREM 9. Let V =3 = {p ∈ C 4 × R 2 | rank(M) = 3} and V 2 = {p ∈ C 4 ×R 2 | rank(M) < 3}, where p = (y 1 , y 2 , z 1 , z 2 , γ 2 , Γ 2 ).
We apply the strategies described in Section 3.

We study the generic case V 2 ∩ V first. This set does cover a dense subset of R 2 . Its intersection with the boundary of B is given by the vanishing of either h 1 or h 2 . The projection on (Γ 2 , γ 2 ) of the set of points of V 2 ∩V such that h 1 = 0 is described by 0 = γ 2 f 2 1 f 2 f 3 which gives us polynomials f 1 , f 2 and f 3 .

The projection on (Γ 2 , γ 2 ) of the set of points of V 2 ∩V such that h 2 = 0 is described by 0 = (2 Γ 2 -γ 2 ) f 1). For both these pairs, there is only 1 singularity in B.

M •   λ 1,1 λ 1,2 λ 2,1 λ 2,2 λ 3,1 λ 3,2 λ 4,1 λ 4,2   = 0 0 0 0 0 0 0 0
with random linear equations ensuring that the matrix (λ i, j ) has rank 2.

Out of the surface γ 2 = 0, this affine variety is a complete intersection (it has dimension 2 and it is given by 9 equations in 11 variables, including the saturation by γ 2 ). The set of critical values of π is described by 0 = (2 Γ 2 -γ 2 ) (Γ 2 + 1) f 2 1 f 2 6 f 2 7 which gives us new polynomials f 6 and f 7 (Γ 2 + 1 has no solutions within our constraint range).

This completes the study of V ∩ V 2 . We now move on to the study of V ∩ V =3 . As described in the algorithm, we define the incidence variety of rank 3 of M, and we saturate successively by the 3-minors of M. Only the first of these subcases is nonempty, and it is described by 0 = (2 Γ 2 -γ 2 ) f 8 f 9 which gives us f 8 and f 9 .

PROOF OF THEOREM 10. Observe first by means of a trivial evaluation that O is a singularity of {D = 0}. We now focus on singularities in B * = B \ {O}. Theorem 9 provides a list of 9 polynomials to which we add our constraints 2 Γ 2 ≥ γ 2 > 0. Let ξ = γ 2 Γ 2 (γ 2 -2 Γ 2 ) ∏ 9 i=1 f i . The complementary of {ξ = 0} is the union of a sequence of connected open semi-algebraic sets where the number of singularities is constant. The routine CylindricalAl-gebraicDecompose of the Maple package RegularChains[Semi-AlgebraicSetTools] provides 1533 sample points. Excluding those at which ξ vanishes and those outside our physical constraints domain, remains a set K c of 570 points. At each point of K c we locate the singularities by computing a Gröbner basis.

We get 187 points of K c such that there exists at least one singularity in B * . We have a set K s of 31 points, each of them corresponding to a couple of ψ-symmetric singularities outside the symmetry plane Π, and a set K p of 156 points corresponding to a unique singularity on Π ∩ B * . For parameters at which ξ does not vanish, the number of singularities in B * is at most two. Points of K s (resp. K p ) are represented in green (resp. blue) in Figs. 1, 3 and 4. Let us evaluate on K c the condition (Γ 2 < 1, f 2 > 0, f 4 < 0) or (Γ 2 > 1, f 2 < 0, f 4 > 0). Indeed the set of points of K c satisfying this condition coincides with K p . This proves item 1). The proof of item 2) is similar.

The general case

The variety V and the semi-algebraic set B are defined as in the previous section. We normalize again by γ 1 = 1, we assume that 2 Γ 1 ≥ 1, 2 Γ 2 ≥ γ 2 > 0, (γ 2 , Γ 2 ) = (1, Γ 1 ), and that Γ 1 = 1, Γ 2 = γ 2 (case of water). THEOREM 11. Splitting the subset of R 3 defined by 2 Γ 2 > γ 2 > 0 and 2 Γ 1 > 1 into open subsets where the number of real singularities of V in the fibers is constant, can be done by cutting out 12 irreducible surfaces, consisting of 5 planes, 3 quadrics, two surfaces of degree 9 and one of degree 14.

These polynomials were obtained by applying the algorithms from Section 3 to our system. The elimination steps were done using both Gröbner bases with FGb or with F5, and with triangular sets with RegularChains. Table 1 presents some timings for these methods (for computations done with interpolation, we give the results as a × b where a is the interpolation degree and b the time taken for each specialized computation). It shows that our algorithm, implemented with either Gröbner bases or regular chains, is more efficient than the direct (general) approach, and it allows to deal with the previously untractable general case. Furthermore, the implementation using Gröbner bases appears to be faster.

DEFINITION 3 .

 3 Let r ∈ {0, . . . , k -1}. The incidence variety of rank r associated with M is the variety V r ⊂ C n+t × (P k-1 (C))

2 Figure 1 :

 21 Figure 1: Curves involved in the definition of the semi-algebraic set G. The blue (resp. green) sample points correspond to pointsin G - 1 ∪ G + 1 (resp. G - 2 ∪ G + 2 ). The circled numbers in each area correspond to the number of singularities in B for parameters in the area. Parameters in the red area are physically irrelevant.

2 Figure 2 :

 22 Figure2: Positions of the parameters corresponding to the pairs water-cerebrospinal fluid (red circle) and water-fat (red square) and the set G (with the same conventions as in Fig.1). For both these pairs, there is only 1 singularity in B.

2 Figure 3 : 2 Γ 2 Figure 4 :

 23224 Figure3: The curves involved in the decomposition of the parameter space (with the same conventions as in Fig.1).

Table 1 :

 1 2 1 f 4 f 5 which gives us new polynomials f 4 and f 5 .Next, we consider the incidence variety V 2 associated with the matrix M: Timings

	System	RegularChain (direct)	Gröbner (direct)	RegularChain (new algo.)	FGb (new algo.)	F5 (new algo.)	CAD
	Water General	1600 s >24 h	100 s >24 h	90 × 200 s	10 s 46 × 200 s	1 s 110 s 4 h (projection step) 50 s

The full results and the source code which produced them are available at http://mercurey.gforge.inria.fr/


A precision analysis of camera distortion models Zhongwei Tang, 1 Rafael Grompone von Gioi, 2 Pascal Monasse, 1 and Jean-Michel Morel 2 1 LIGM, UMR 8049, École des Ponts, UPE, Champs-sur-Marne, France 2 CMLA, ENS-Cachan, Cachan, France

Abstract-This paper addresses the question of identifying the right camera direct or inverse distortion model permitting a high subpixel precision fit to real camera distortion. Five classic camera distortion models are reviewed and their precision compared for direct or inverse distortion. By definition, the three radially symmetric models can only model a distortion radially symmetric around some distortion center. They can be extended to deal with non-radially symmetric distortions by adding tangential distortion components, but still may be too simple for very accurate modeling of real cameras. The polynomial and the rational models instead miss a physical or optical interpretation, but can cope equally with radially and nonradially symmetric distortions. Indeed, they do not require the evaluation of a distortion center. When requiring high precisions, we found that the distortion modeling must also be evaluated primarily as a numerical problem. Indeed, all models except the polynomial involve a non-linear minimization which increases the numerical risk. The estimation of a polynomial distortion model leads instead to a linear problem, which is secure and much faster. We concluded by extensive numerical experiments that, although high degree polynomials were required to reach a high precision of 1/100 pixels, such polynomials were easily estimated and produced a precise distortion modeling without over-fitting. Our conclusion is validated by three independent experimental setups: The models were compared first on the lens distortion database of the Lensfun library by their distortion simulation and inversion power; second by fitting real camera distortions estimated by a non parametric algorithm; and finally by the absolute correction measurement provided by photographs of tightly stretched strings, warranting a high straightness.

Index Terms-distortion measurement, camera calibration

I. INTRODUCTION

The pinhole camera model is widely used in computer vision applications because of its simplicity and its linearity in terms of projective geometry [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. But real cameras deviate from the ideal pinhole model. The main geometric deviation is a lens geometric distortion [START_REF] Brown | Close-range camera calibration[END_REF], possibly complicated by a deviation from planarity of the CCD shape. Thus an accurate camera distortion correction is the first step towards high precision 3D metric reconstruction from photographs. With the steady progress in lens quality and camera resolution, highprecision 3D reconstructions become feasible. But they require in turn higher camera distortion precisions than those provided by classic methods. The object of this paper is to investigate the validity of distortion models at the light of precision requirements increased by two or three orders of magnitude. This increased accuracy requires a new methodology for evaluating distortion models. Five models are studied and evaluated in the paper: the radial [START_REF] Brown | Close-range camera calibration[END_REF], division [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF], FOV [START_REF] Devernay | Straight lines have to be straight[END_REF], polynomial (such as bicubic [START_REF] Kilpelä | Compensation of systematic errors of image and model coordinates[END_REF]), and rational [START_REF] Claus | A rational function lens distortion model for general cameras[END_REF], [START_REF] Hartley | The cubic rational polynomial camera model[END_REF] models. They rely on several different hypotheses on the underlying distortion model. Clearly their precision depends on which model hypothesis is valid. Radially symmetric models are very precise for radially symmetric distortions. They can be extended to treat non-radially symmetric distortions by adding tangential components. The polynomial model and the rational model impose the fewer constraints on the distortion, to the cost of an increased number of model parameters. Thay way they can cope with radial and non-radial symmetric distortions as well. Among them the polynomial is the only linear model. Its estimation is a simple matrix inversion. All the other models require a somewhat complex non-linear minimization. We shall see that the linearity of the polynomial model reduces its numerical risk and makes its modeling precision closer to 1/100 pixels for realistic camera distortions.

It could be argued that a correct model should be based on physical measurements on systems of lenses. Surprisingly enough, there is little physical background for the distortion models in the literature. According to Weng et al. [START_REF] Weng | Camera calibration with distortion models and accuracy evaluation[END_REF], lens distortion can be decomposed into three effects: radial distortion, decentering distortion and thin prism distortion. Nevertheless, it is only marginally based on a physical background, and it is not clear that the real distortion precisely follows this model. In fact, the final distortion is the result of the cumulated effects of a complex lens system, of the camera geometry, and of the (not perfectly planar) shape of the image sensor. One is therefore led to figure out a flexible model with enough parameters to simulate any plausible distortion. In absence of a physical model, the model classification approach adopted here will be to look for models that actually cope with any realistic distortion, at a given precision.

Many works in the computer vision community assume that the distortion is radially symmetric around the center of the image, and the radial, division and FOV models are often used to simulate and inverse the radially symmetric distortion. Even though these models are not exactly (algebraically) invertible, the simulation and correction can be obtained at very high precision if the order of the model is high enough, as we will show in Section III. This explains why these models in the literature are used with the interchangeable roles of distorted points and undistorted points. For example, direct distortion models are used in global camera calibration [START_REF] Tsai | A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses[END_REF], [START_REF] Zhang | A flexible new technique for camera calibration[END_REF], [START_REF] Lavest | Do we really need accurate calibration pattern to achieve a reliable camera calibration[END_REF], [START_REF] Weng | Camera calibration with distortion models and accuracy evaluation[END_REF]. Yet, in most plumb-line methods [START_REF] Brown | Close-range camera calibration[END_REF], [START_REF] Devernay | Straight lines have to be straight[END_REF], [START_REF] Alvarez | An algebraic approach to lens distortion by line rectification[END_REF], [START_REF] Alvarez | Algebraic lens distortion model estimation[END_REF], [START_REF] Santana-Cedrés | An iterative optimization algorithm for lens distortion correction using two-parameter models[END_REF], [START_REF] Prescott | Line-based correction of radial lens distortion[END_REF], [START_REF] Pajdla | Correcting radial lens distortion without knowledge of 3-D structure[END_REF], [START_REF] Claus | A plumbline constraint for the rational function lens distortion model[END_REF] or some pattern-free methods [START_REF] Stein | Lens distortion calibration using point correspondences[END_REF], [START_REF] Zhang | On the epipolar geometry between two images with lens distortion[END_REF], [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF], [START_REF] Li | A non-iterative method for correcting lens distortion from nine point correspondences[END_REF], [START_REF] Thirthala | Multi-view geometry of 1d radial cameras and its application to omnidirectional camera calibration[END_REF], [START_REF] Claus | A rational function lens distortion model for general cameras[END_REF], [START_REF] Pajdla | Automatic generator of minimal problem solvers[END_REF], [START_REF] Byrod | Fast and robust numerical solutions to minimal problems for cameras with radial distortion[END_REF], [START_REF] Kukelova | A minimal solution to the autocalibration of radial distortion[END_REF], the very same radial correction models are used without any fuss to inverse the distortion.

In practice, the distortion center is unknown and we do not even a priori know if the distortion is radially symmetric around a certain center. So the radially symmetric models with an arbitrarily fixed distortion center (typically, the image center) is not always a good model. The common practice is to fix the distortion center at some reasonable position and add tangential distortion components to the model. All the global camera calibration methods adopt this strategy by fixing the distortion center at the principal point. An alternative is to estimate the distortion center that makes the radially symmetric models fit optimally to the distortion. For example, Hartley and Kang [START_REF] Hartley | Parameter-free radial distortion correction with centre of distortion estimation[END_REF] propose to estimate the distortion center under the assumption that the distortion is radially symmetric and monotone. Being more general, the polynomial model and the rational model are invariant to the distortion center and can be directly used to model the distortion without estimating any distortion center.

Our aim is to find models which are flexible and easy to use, producing small residual error no matter what reasonable distortion has been applied. Ideally we would like to test on all possible distortions of any existing camera lens. Since exhaustive testing is impossible, we resorted to the lens database of the LensFun library, which is known to contain the most complete freely available lens profiles. All the distortions provided by Lensfun are radially symmetric around the image center.

The question we raise here only makes sense within fixed accuracy bounds. This is a first caveat: any evaluation must be performed at a given precision. As a matter of fact, for off-the-shelf cameras, most distortion models perform well at a 1 pixel precision. The question is different when we aim at sub-pixel precisions. These precisions, up to 1/100 pixels, are highly desirable when using cameras for stereovision or photogrammetric tasks. Indeed, in stereovision a determining factor is b/h, the ratio of baseline to average depth. A high ratio improves accuracy of depth determination by triangulation, but makes automatic image matching more difficult because of the significant viewangle change. On the contrary, a low ratio eases the point-matching step but requires sub-pixel estimation to avoid a strong quantization of recovered depths. The latter case is advisable for automatic processing, hence the need for highly accurate distortion correction.

The other caveat is that, although distortion models reflect a model of the optical lens, the real distortion must actually involve the whole system lens plus CCD. There is no way to guarantee that a CCD is absolutely flat, or exactly perpendicular to the optical axis. This explains why the camera distortion modeling remains, after all, an empirical question where no physical argument can be final. The ultimate decision is numerical.

The various distortion models will be carefully compared on lens distortions in the Lensfun library, permitting to quantify the ideal attainable precision. Then, the same models will be compared on their capacity to fit to real camera distortions (estimated by a non-parametric algorithm [START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF]). Finally, the distortion correction accuracy by each model will be evaluated by using the plumb-line approach, with photographs of tightly stretched strings, warranting a high straightness, and giving absolute measurements of the correction quality [START_REF] Tang | Highprecision camera distortion measurements with a "calibration harp[END_REF]. In short, there will be three different numerical validations of our conclusions.

This paper is organized as follows. Section II reviews five classic distortion models. Their power to simulate and inverse distortions are evaluated in Section III by synthetic experiments. Section IV and V describe the experiments performed on real cameras. Section VI sums up the lessons learned from these experiments.

II. DISTORTION AND CORRECTION MODELS

We denote by (x u , y u ) the coordinates of an undistorted point as would be observed by an ideal pinhole camera. Due to the lens geometric distortion, this point will be observed at coordinates (x d , y d ). The distortion is modeled by a function f that transforms undistorted to distorted coordinates,

x d = f x (x u , y u ) y d = f y (x u , y u ). (1)
A correction model g performs the inverse transformation,

x u = g x (x d , y d ) y u = g y (x d , y d ). (2)
A particularly interesting case is when the function f or g shows radial symmetry relative to a fixed distortion center (x c , y c ). In that case we obtain a compact formulation using radial coordinates xu = x u -x c , ȳu = y u -y c , xd = x d -x c and ȳd = y d -y c ; then, the distortion can be expressed as the transformation of the undistorted radius r u = x2 u + ȳ2 u to the distorted radius r d = x2

d + ȳ2 d . We start by reviewing the most current models. To treat both directions in a neutral way, we will write the models as transforming from coordinates (x 1 , y 1 ) to (x 2 , y 2 ). When a model would be used as a distortion model, (x 1 , y 1 ) will correspond to (x u , y u ) and (x 2 , y 2 ) to (x d , y d ), and it is the opposite when used as a correction model.

The radial model displaces a point along its radial direction originating at the distortion center. The new radius r 2 is a function of the original radius r 1 ,

r 2 = r 1 (k 0 + k 1 r 1 + k 2 r 2 1 + • • • ). (3)
The parameter k 0 represents a scaling and therefore does not introduce distortion. The scaled image is distorted by k 1 , k 2 , . . . If all are positive, we have a pincushion distortion; if all are negative, a barrel distortion. Mustache distortion occurs if the signs are not the same 1 . Note that the distortion center (x c , y c ) is also a parameter of radial models. The division model is obtained by simply putting the factor term

(k 0 + k 1 r 1 + k 2 r 2 1 + • • • ) on the denominator, r 2 = r 1 k 0 + k 1 r 1 + k 2 r 2 1 + • • • . ( 4
)
In these models, high-order coefficients are needed to model extreme distortion in fish-eye lenses or other wide angle lens systems. A more sparse representation is obtained by parameterizing the distortion by the field of view (FOV).

The only parameter of the FOV model is the field of view parameter ω:

r 2 = tan(r 1 tan(ω)) tan(ω) . 2 (5)
The Taylor expansion of the FOV model around r 1 = 0 is

r 2 = r 1 + tan 2 (ω) 3 r 3 1 + 2 tan 4 (ω) 15 r 5 1 + • • • . (6)
So the FOV model is a radial model with odd terms around r 1 = 0. Note that the order-0 scale coefficient is fixed to be 1 and all the other coefficients are coupled to ω. So to model more complex distortions, the authors proposed to complete it with k 0 , k

1 , k 3 , • • • : 3 r 2 = tan(r 1 tan(ω)) tan(ω) + r 1 (k 0 + k 1 r 1 + k 3 r 3 1 • • • ). (7)
In the polynomial model the distortion is modeled as a polynomial in x 1 and y 1 . For example, the third order (bicubic) polynomial model is

x 2 = a 1 x 3 1 + a 2 x 2 1 y 1 + a 3 x 1 y 2 1 + a 4 y 3 1 + a 5 x 2 1 + a 6 x 1 y 1 + a 7 y 2 1 + a 8 x 1 + a 9 y 1 + a 10 , y 2 = b 1 x 3 1 + b 2 x 2 1 y 1 + b 3 x 1 y 2 1 + b 4 y 3 1 + b 5 x 2 1 + b 6 x 1 y 1 + b 7 y 2 1 + b 8 x 1 + b 9 y 1 + b 10 , (8)
depending on the (3 + 1)

• (3 + 2) = 20 parameters a 1 , • • • , a 10 , b 1 , • • • , b 10 .
More generally, a polynomial model of order n depends on (n+1)•(n+2) parameters. The rational function model is a quotient of two polynomials. A second order rational function model can be written as

x 2 = a 1 x 2 1 + a 2 x 1 y 1 + • • • + a 5 y 1 + a 6 c 1 x 2 1 + c 2 x 1 y 1 + • • • + c 5 y 1 + c 6 , y 2 = b 1 x 2 1 + b 2 x 1 y 1 + • • • + b 5 y 1 + b 6 c 1 x 2 1 + c 2 x 1 y 1 + • • • + c 5 y 1 + c 6 . (9)

III. PRECISION EVALUATION

We shall evaluate the distortion simulation and inversion power of the above mentioned models. Being a theoretical property of model families, the precision in both directions can be genuinely evaluated by synthetic experiments.

The test begins by synthesizing realistic distortions produced by a camera lens. Since it is impossible to exhaustively obtain all existing lens distortion profiles, we resort to the lens database of the public library LensFun, which inherits the database of the commercial software PTLens and has the most comprehensive freely available lens database (see Table I for a few lens examples. A complete lens list supported by Lensfun can be found at http://lensfun.sourceforge.net/lenslist/. We use version 0.3.2 of LensFun comprising more than 3500 models.). In Lensfun, the distortion is calibrated with some predefined models (see Table II), based on the matching points between two images taken by the same camera on the same focal length. 4 The final calibrated distortion models in Lensfun are represented in the normalized image domain II.

[-1.0, +1.0] × [-1.0, +1.0] 5 .
Model Formulation ptlens

r d = ru(1 -a -b -c + cru + br 2 u + ar 3 u ) poly3 r d = ru(1 -k 1 + k 1 r 2 u ) poly5 r d = ru(1 + k 1 r 2 u + k 2 r 4 u )
TABLE II: Models used to generate distortion in Lensfun.

According to each calibrated lens model in LensFun, we can synthetically generate distorted/undistorted point pairs in the normalized image domain. The simulation and correction precisions are then verified by identifying the best parameters through (1) and ( 2) respectively. In other words, both (x u , y u ) and (x d , y d ) are known in the synthetic tests, and the question is how well the distorted points (x d , y d ) can be approached by f x (x u , y u ) and f y (x u , y u ), and how well the ideal points (x u , y u ) can be approached by g x (x d , y d ) and g y (x d , y d ). We want to compute the coefficients of f x and f y by minimizing the difference between the observed distortion and the simulated distortion. The energy to be minimized can be written as

C a = f x (x u , y u ) -x d 2 + f y (x u , y u ) -y d 2 dx u dy u .
(10) Similarily, the energy to be minimized for estimating the coefficients of g x and g y can be written as

C c = g x (x d , y d ) -x u 2 + g y (x d , y d ) -y u 2 dx d dy d . (11)
In practice, we generated a total of 2M pairs of distorted/undistorted points {(x di , y di ), (x ui , y ui )} i=1,...,2M uniformly distributed in the normalized image domain. Among these pairs, M of them were used to estimate the parameter by minimizing the discrete energy

D a = M i=1 f x (x ui , y ui ) -x di 2 + f y (x ui , y ui ) -y di 2 (12)
and

D c = M i=1 g x (x di , y di ) -x ui 2 + g y (x di , y di ) -y ui 2 (13)
respectively. In order to verify the precision of the tested models, the estimated parameters were then used to compute the average residual error: The unknown parameters include the specific parameters of each distortion model and the distortion center (x c , y c ). This center (x c , y c ) is irrelevant for the polynomial model and the rational function model, both models being invariant to the distortion center. Yet for the other models, (x c , y c ) is an important parameter that must be accurately estimated to achieve precise distortion modeling [START_REF] Clarke | The principal point and CCD cameras[END_REF], [START_REF] Hartley | Parameter-free radial distortion correction with centre of distortion estimation[END_REF].

Linear or non-linear optimization methods are required to minimize the energy in ( 12) and [START_REF] Hartley | The cubic rational polynomial camera model[END_REF]. For the radial model (if the distortion center (x c , y c ) is known) and the polynomial model, the energies can be minimized by solving a linear system Ak = b,

where A is a matrix containing the different powers for all the distorted (or undistorted) points, k is formed with the correction (or distortion) model coefficients, and b contains the undistorted (or distorted) coordinates. The set of model coefficients with least error is obtained by minimizing the residual Akb 2 , which results in

k = A T A -1 A T b. ( 15
)
Since the synthetic tests can be done in the normalized image domain, the coefficient matrix A is usually well conditioned. For all other models, a non-linear method must be used, even if (x c , y c ) is known. The minimization was performed by an incremental Levenberg-Marquardt (LM) algorithm which estimates the parameters in increasing order: The algorithm starts estimating the parameters of a second-order model; the result is used to initialize the model with the next higher order, and the process continues until the aimed order. The Jacobian matrix J required by LM was computed explicitly to make the algorithm efficient and more precise than with finite differences. We chose quite strict LM stopping criteria: each individual LM stops after 1000 iterations or when the energy does not decrease by more than 1e-10. Even though this strategy was complex, it avoided some local minima and is safer than performing LM directly on the model at the aimed order.

Parameters initialization: For non-linear minimization, it is important to initialize the parameters with reasonable values. The common practice for the radial model and the division model is to initialize the 0th-order parameter k 0 to be 1 and all the other parameters to be 0. The parameter ω in the FOV model is interpreted as half the field of view of the camera. For ordinary cameras, a typical value of ω is about π 6 ≈ 0.52, which can be can be used as initial value.

The initialization of the rational function model is more involved. As shown in [START_REF] Claus | A rational function lens distortion model for general cameras[END_REF], it can be solved linearly by using a "lifted process" technique. More precisely, (x 1 , y 1 ) and (x 2 , y 2 ) are related by a 3 × 6 matrix A:

λ(x 2 , y 2 , 1) T = A(x 2 1 , x 1 y 1 , y 2 1 , x 1 , y 1 , 1) T , (16)
A =   a 1 a 2 a 3 a 4 a 5 a 6 b 1 b 2 b 3 b 4 b 5 b 6 c 1 c 2 c 3 c 4 c 5 c 6   . (17)
A can be solved linearly by using at least 9 distorted/undistorted point pairs, similar to the classic homography estimation. This linear solution is used in the experiments to initialize the incremental LM algorithm, since it minimizes the algebraic error, which is not directly related to the geometric error we want to minimize and sometimes leads to undesirable result.

Distortion center initialization:

The distortion center is a sensitive parameter for distortion models and should be accurately estimated to model or remove distortion precisely [START_REF] Clarke | The principal point and CCD cameras[END_REF]. The common practice is to put the initial distortion center at the center of image. However, this is not always a safe initialization because the true solution can be away from the image center up to 10% of the image domain [START_REF] Clarke | The principal point and CCD cameras[END_REF]. A more delicate approach was proposed [START_REF] Hartley | Parameter-free radial distortion correction with centre of distortion estimation[END_REF] to linearly estimate the distortion center for radially symmetric models, which can be used as an initialization for our non-linear minimization. Based on the monotonicity assumption of the radially symmetric distortion, this method exploits the observation that the distorted points are obtained by moving the undistorted points away from (or to) the distortion center, similar to the motion of points seen by a camera moving forward toward/backward from a scene. So the distortion center plays the same role as the epipole of the fundamental matrix in a special situation, where the camera undergoes a pure translation between two views [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. This observation leads to the following epipolar geometry [START_REF] Hartley | Parameter-free radial distortion correction with centre of distortion estimation[END_REF]:

x T di [e] × x ui = x T di ([e] × H)x ci = x T di Fx ci = 0, (18)
with [e] × the skew-symmetric 3 × 3 matrix representing the cross product by the distortion center e. The undistorted points x ui are projected from the 3D points x ci on a planar calibration pattern through a pinhole camera:

x ui = Hx ci .
The matrix F can be linearly computed from at least 7 pairs of {x di , x ui } correspondences [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. In our case, H is the 3×3 identity matrix and we have simply: x T di [e] × x ui = 0. The distortion center e can then be computed as the right null vector of the matrix [e] × , which can be solved by at least 7 distorted/undistorted point pairs, just like the fundamental matrix estimation in binocular stereo [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

A. Experiments with known distortion center

We first consider the case where the distortion center is known. Then the radial model and the polynomial can be estimated by a matrix inversion, whereas all the other models require an incremental LM minimization.

Simulation: Notice that the distortions are all generated by radially symmetric models, see Table II. So the radial model should be able to ideally simulate these distortions if it contains all the terms used in the distortion generation in Table II. Fig. 1a shows that the simulation precision is close to the machine limit except for the curve corresponding to the third-order radial model, which cannot simulate the distortions generated by higher order models (Table II). The division model can also simulate the distortion well. Even though this simulation is not exact, the precision increases with the order of the division model. In fact, a radially symmetric distortion can be written as

r 2 = r 1 (k 0 + k 1 r 1 + k 2 r 2 1 + • • • ) = r 1 1 (k0+k1r1+k2r 2 1 +••• ) = r 1 1 k0 -k1 k 2 0 r 1 + 1 k0 ( -k2 k0 + k 2 1 k 2 0 )r 2 1 + • • • , ( 19
)
which means that the radially symmetric distortion can be well simulated by the division model with high enough order.

For the FOV model, as we have shown before, since it has a high correlation with second order coefficient k 2 and has been completed by the other coefficients, it should also be able to simulate the radially symmetric distortion very precisely. However, in practice, both the division model and the FOV model estimation require a non-linear minimization and the resulting numerical risk explains why we obtained a less precise simulation in practice (Fig. 1b and1c).

As for the polynomial model, it can simulate ideally the radially symmetric distortions composed of only even-order terms

r 2 = r 1 (k 0 + k 2 r 2 1 + k 4 r 4 1 + • • • ) because x 2 = x 1 (k 0 + k 2 r 2 1 + • • • ) = k 0 x 1 + k 2 x 3 1 + k 2 x 1 y 2 1 + • • • y 2 = y 1 (k 0 + k 2 r 2 1 + • • • ) = k 0 y 1 + k 2 x 2 1 y 1 + k 2 y 3 1 + • • • are just
two polynomials, with some monomials missing though. If odd-order terms also appear:

r 2 = r 1 (k 0 + k 1 r 1 + k 2 r 2 1 + • • •
), the polynomial model cannot ideally simulate the distortion. But the simulation precision increases with the order of polynomial model. Two polynomial models of adjacent orders (for example, the third-order and fourth order polynomials) have very similar performance due to the fact that the distorted/undistorted point pairs are almost radially symmetric and thus the odd-order terms in the polynomial model play a more important role to model such points. Also the error curve of the polynomial model is cut into three pieces due to the fact that the synthesized distortions can be classified into three groups: one group with only odd-order terms, one group with odd-order terms and one even-order term, and one group with odd-order terms and two even-order terms. These observations explain what we observe in Fig. 1e. In fact, if the radial model contains only odd terms, we observe the similar curves in Fig. 1d as in Fig. 1e. The rational function model is more general than the polynomial model. So it should in theory improve the simulation power of the polynomial model (Fig. 1f).

The maximum residual error among the points, represented as dashed lines in the figure, is typically higher by a factor 10 than the average error. This shows that even the worst error remains low for most camera lenses.

Inversion: The correction power is the precision achieved by the distortion models to inverse the distortions in Lensfun. For any tested model, the inversion is not exact for the radially symmetric distortions in Lensfun. However, to inverse the distortion following

r 2 = r 1 (k 0 + k 1 r 1 + k 2 r 2 2 + • • • ), we have r 1 = r 2 k 0 + k 1 r 1 + k 2 r 2 1 + • • • = r 2 1 k 0 - k 1 k 2 0 r 1 + • • • = r 2 1 k 0 - k 1 k 2 0 • r 2 k 0 + k 1 r 1 + k 2 r 2 1 + • • • + • • • = r 2 1 k 0 - k 1 k 2 0 • r 2 ( 1 k 0 - k 1 k 2 0 r 1 + • • • ) + • • • = r 2 1 k 0 - k 1 k 3 0 r 2 + k 2 1 k 5 0 r 2 2 + • • • (20)
or

r 1 = r 2 k 0 + k 1 r 1 + k 2 r 2 1 + • • • = r 2 k 0 + k 1 • r2 k0+k1r1+k2r 2 1 +••• + • • • = r 2 k 0 + k 1 r 2 1 k0 -k1 k 2 0 r 1 + • • • + • • • = r 2 k 0 + k1 k0 r 2 - k 2 1 k 3 0 r 2 2 • • • , (21)
which means that radially symmetric distortions are invertible by the radial model, the division model and the FOV model. The inversion is not exact and the precision increases with the order of the inversion model. In addition, with the same argument as before, the polynomial model and the rational function model can also be inverse of the radially symmetric distortion, in particular when the radially symmetric distortion is composed of only even-order terms. All the curves are shown in Fig. 1g-1l. It was shown in [START_REF] Santana-Cedrés | An iterative optimization algorithm for lens distortion correction using two-parameter models[END_REF] that the second order division model is superior to the radial of the same order; our experiments do not show such advantage for third or higher orders. Today's digital cameras produce images of millions of pixels with a minimal order of magnitude of 1000 × 1000. We aim at an average precision about 10 -2 pixel for images of such size. This is about 10 -5 in the normalized image domain. All the models are able to simulate and correct the radially symmetric distortions provided by the Lensfun library at this precision if the order is high enough. In this restricted setting, they are all equivalent, notwithstanding the difficulty of minimization for some.

B. Experiments with unknown distortion center

In practice, the distortion center (x c , y c ) is unknown and we have to estimate it as well. The distorted/undistorted point pairs were generated in the same way as before except that the distortion center (x c , y c ) was unknown. We used the linear method proposed in [START_REF] Hartley | Parameter-free radial distortion correction with centre of distortion estimation[END_REF] to first obtain an initial estimate of the distortion center, which was then refined with the (l) rational model for correction Fig. 1: The average residual error (dahed line: maximal error) of different models when used to simulate the distortion or its inverse (i.e., the correction) synthesized following the models and the parameters provided by the Lensfun library. The horizontal axis is the lens index (or one lens at different focal lengths) and the vertical axis is the average residual error represented in log 10 in the normalized domain. The average residual errors are ordered increasingly. This figure is better viewed electronically in color. other parameters of distortion models in the incremental LM minimization. Since the initial estimate of (x c , y c ) is already precise, we finally obtained results very similar to the models with (x c , y c ) known and fixed. For the polynomial model and the rational model, since they are invariant to the distortion center, it is not suprising that they fall into the same precision. We do not show these very similar curves to make the paper more concise.

C. Comparison of radial and division models

Based on observations of Fitzgibbon [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF], the division model is assumed to be more precise than the radial model, at the same number of parameters. This is confirmed by the experiments of Santana-Cedrés et al. [START_REF] Santana-Cedrés | An iterative optimization algorithm for lens distortion correction using two-parameter models[END_REF]. The compared models are infinitely differentiable, with one or two parameters, that is (3) or ( 4) with: k 0 = 1, k 1 = k 3 = 0, and k 2 , possibly k 4 , as variables. To test them on LensFun models, we need to let k 0 as variable, as LensFun models are not necessarily close to the identity near the origin. In [START_REF] Santana-Cedrés | An iterative optimization algorithm for lens distortion correction using two-parameter models[END_REF], such radial and division models are used to correct the distortion. Therefore, in Figure 2, we measure the precision obtained when trying to inverse the distortions of LensFun models. The division model is indeed more precise than the radial model with the same number of parameters, yet the improvement is marginal. These experiments do not support the observations of the cited articles.

D. Realistic distortion

We remark that the lens profile in the LensFun database is only approximatively calibrated: The provided lens distortions are all radially symmetric around the image center. This is the drawback of most commercial software and explains why they correct the distortion at a precision no better than 0.5 pixel. In reality, we are not sure if the real distortion is symmetric around a certain distortion center. The model complemented by the tangential distortion components is perhaps closer to the real distortion [START_REF] Weng | Camera calibration with distortion models and accuracy evaluation[END_REF]:

xu = xd + xd (k 1 r 2 d + • • • + k 5 r 10 d ) + p 1 r 2 d + 2x 2 d + 2p 2 xd ȳd , ȳu = ȳd + ȳd (k 1 r 2 d + • • • + k 5 r 10 d ) + p 2 r 2 d + 2ȳ 2 d + 2p 1 xd ȳd , (22)
with p 1 , p 2 , p 3 parameters for decentering distortion and s 1 , s 2 parameters for thin prism distortion. They contribute to both radially symmetric distortion and tangential distortion. This model is often adopted in bundle adjustment algorithms. The parameters calibrated with Lavest et al. algorithm [START_REF] Lavest | Do we really need accurate calibration pattern to achieve a reliable camera calibration[END_REF] on an EF-S 18-55mm f/3.5-5.6 lens mounted on a Canon EOS 30D camera were used to synthesize the distortion according to [START_REF] Pajdla | Correcting radial lens distortion without knowledge of 3-D structure[END_REF]. The synthesized distortion was then corrected by the different models through incremental LM minimization.

The results are shown in Table III. The radial+tangential model has the distortion center fixed at the known distortion center and inverses the distortion synthesized by the same model. The parameters p 1 and p 2 were initialized at 0. For three radially symmetric models, we estimated also the distortion center which makes them best fit the distortion. The polynomial model and the rational model are invariant to the distortion center. It seems that the radial+tangential model can better inverse the distortion (fitting error 0.005 pixels) than the three radially symmetric models (fitting error 0.01 pixels). In fact, since the synthesized distortion follows the above model in [START_REF] Pajdla | Correcting radial lens distortion without knowledge of 3-D structure[END_REF], it can be inversed by a model of the same family if the tangential component and the radial component satisfy the following constraints (see Fig. 3a):

∆t = ∆d∆t ∆u , ∆r = ∆r 2 + ∆t 2 -∆t 2 .
∆d is the distance between the distorted point to the distortion center and ∆u is the distance between the undistorted point to the distortion center. ∆r and ∆t are the radial component and the tangential component of the synthesized distortion. ∆r and ∆t are the radial component and the tangential component of inversed distortion. Similarly, the synthesized distortion can also be modeled by the radially symmetric models if the following constraint is satisfied (see Fig. 3b):

∆r = ∆r 2 + ∆t 2 (23)
where (x c , y c ) denotes the new distortion center for the radially symmetric models. These constraints are complex and no closed-form formula exists. They can be simulated but not exactly satisfied. This explains why the polynomial and the rational models are even more precise and yield errors on the order of 1e-5 pixels. But we should remember that the distortion is synthesized following [START_REF] Pajdla | Correcting radial lens distortion without knowledge of 3-D structure[END_REF] and we are not sure if the real distortion can be precisely modeled in this way. We shall confirm the precision and the stability of the models under consideration in the next real experiments. We will see in the real experiments that the radial+tangential models perform sometimes worse than the radially symmetric models. This poses the problem of model selection in practice. but could be performed on any distortion model obtained by blind correction. This method requires a highly textured planar pattern, which is obtained by printing a textured image and pasting it on a very flat object (a thick, heavy and rigid aluminium plate was used in the experiments). Another option is to use a high density grid pattern as in [START_REF] Lee | Accurate determination of distortion for smartphone cameras[END_REF]. Two photos of the pattern were taken by a same camera. The distortion is estimated (up to a homography) as the diffeomorphism mapping the original digital pattern to its photograph. The algorithm is summarized below.

1) Take two slightly different photographs of a textured planar pattern (Fig. 4) with a camera whose settings are manually fixed (disable automatic mode); 2) apply the SIFT method [START_REF] David G Lowe | Distinctive image features from scale-invariant keypoints[END_REF] between the original digital pattern and two photographs respectively, to obtain matching pairs of points; 3) triangulate and interpolate the SIFT matches from the digital image to two photographs respectively; 4) use a loop validation to eliminate the outlier matches from the digital pattern to one of the two photographs; 5) use a vector filter to remove the few remaining outliers matches from the digital pattern to that photograph; 6) refine the precision of the SIFT matching by correcting each matching point in one image with the local homography estimated from its neighboring matching points; 7) triangulate and interpolate the refined inlier matches to get a dense reverse distortion field from the digital pattern to that photograph; 8) apply the reverse distortion field to any image produced by the same camera to correct the distortion. For more details please refer to [START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF]. The matching pairs delivered by step 6 (about 8000 in our experiments) in the algorithm above do not contain gross outliers and are precise thanks to the local filtering. So we could directly try all the models to fit these matchings. The residual fitting error shows to what extent the models are faithful to a real camera distortion. Under the assumption that the textured pattern is flat, the mapping from the digital pattern to the photo can be modeled as SDH, with H the homography from the digital pattern to the photo, D the non-linear lens distortion and S a diagonal matrix to model the slant of the CCD matrix:

H =   h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 1   , S =   1 0 0 0 s 0 0 0 1   . ( 24
)
Since the polynomial and the rational function models can simulate well H and S, we can use these models to simulate the distortion field without explicitly estimating the homography. Nevertheless, for the radially symmetric models, it is necessary to take into account H and S when approximating the distortion. Indeed H and S are generally not radially symmetric. Thus, we have 9 additional parameters to estimate, besides the parameters of radially symmetric distortion models and their distortion center. The radial+tangential model is also tried with H and S (with the distortion center fixed at the center of image). The polynomial model can again be solved linearly. For all the other models, an incremental LM minimization was used to estimate the distortion center, the distortion parameters, H and S. The matrix H was initialized as the homography linearly estimated from the digital pattern to the photo and s is initialized at 1. The other parameters were initialized in the same way as we did for the synthetic tests. We worked in the normalized image domain to avoid possible numerical problems.

Half the matching pairs were used to estimate the parameters for the different models, and the other half to evaluate the average fitting error. We tried two cameras: a Canon EOS 30D with EF-S 18-55mm lens and a Canon EOS 40D with EF 24-70mm lens. Both extreme focal lengths were tested: 18mm and 55mm for EF-S 18-55mm lens, 24mm and 70mm for EF 24-70mm lens. The results are recapitulated in Table IV. They show that by combining H and S to model the inclination between the camera and the pattern, all of the radially symmetric models give almost the same fitting error, which becomes stable with the increase of the model order. The similar estimation of H and S indicates that the minimization process is stable. The radial+tangential model is not always better than the radially symmetric models.

Sometimes, it is even slightly worse (24mm focal length, Table IV).

The polynomial and the rational function models give a stable fitting error, which is smaller than the radially symmetric models. The stability of the fitting error confirms that none of the models suffers from numerical instability or noise fitting. We remark that the residual with the rational function model becomes stable just a little faster than the polynomial model. However, solving the rational function model requires a complex incremental LM minimization, which is a timeconsuming process requiring a good initialization and does not always ensure the convergence to the global minimum, whereas the polynomial model can always be solved linearly.

V. PLUMB-LINE VALIDATION

It should be noted that the non-parametric method does not give a ground truth. It is just a non-parametric estimation of the camera distortion, and it is subject to errors. Thus, we needed a more objective evaluation to check the quality of the correction models. To this purpose, a physical frame with tightly stretched opaque fishing strings was built [START_REF] Tang | Highprecision camera distortion measurements with a "calibration harp[END_REF]. The physical tension of the strings guarantees a very high straightness. Once the parameters of the models are estimated, a distortion field can be constructed and applied for the distortion correction of images of strings taken by the same camera with the same fixed lens configuration, see Fig. 5. The "average straightness error" is the average distance from the sub-pixel edge points (computed by an improved Canny filter [START_REF] Devernay | A non-maxima suppression method for edge detection with sub-pixel accuracy[END_REF]) of the corrected lines to the corresponding regression line (see [START_REF] Tang | Highprecision camera distortion measurements with a "calibration harp[END_REF] for more details). Table V recapitulates the average distance measured on the strings image corrected by different models for the lens we tested. The stabilized straightness error implies again that none of the models has the problem of noise fitting, which guarantees the correction quality and stability. The radial+tangential model does not always perform better than the radially symmetric models and sometimes it is even worse (focal length 24mm). The radially symmetric models with varying distortion center seems enough to fit the real camera distortion precisely when the distortion is not very big. For the big distortion with short focal length 18mm, the polynomial model and the rational model are twice more precise.

We also remark that the measurement precision provided by the fishing strings is limited by the detection precision of edge points and the small intrinsic oscillation of the fishing strings. The limit measurement precision is around 0.03 ∼ 0.04 pixels, which appears many times in Table V. This means that the distortion was corrected at approximately the limit precision that could be measured. This fact might also explain why the polynomial and the rational models giving smaller fitting error do not always lead to smaller straightness error. Some other more precise instruments might be able to measure the different errors obtained by the two models.

VI. CONCLUSION

We evaluated and compared the precision of five distortion models. The radially symmetric models are ideal for modeling the radially symmetric distortion, but in practice, the related non-linear optimization raises a numerical risk. These models can be extended by adding tangential distortion components. However, this is not always a good choice and can pose the problem of model selection. The polynomial and the rational models are more universal in the sense that they are invariant to the distortion center and can model realistic non-radially symmetric distortion. The polynomial model requires a little higher degrees than the rational model to achieve the precision on the order of 1/100 pixels. This high degree raises no [START_REF] Rudakova | Camera matrix calibration using circular control points and separate correction of the geometric distortion field[END_REF]. In fact, the results of all the tested parametric models are visually identical, according to the accuracy shown in Tables IV andV.

computational problem for the polynomial model, which can be solved linearly after a correct conditioning.

The polynomial and the rational models are not adapted to global bundle adjustment, where the internal/external parameters and the distortion model are estimated simultaneously. The distortion correction must be dealt with as an independent and previous step to camera calibration [START_REF] Rudakova | Camera matrix calibration using circular control points and separate correction of the geometric distortion field[END_REF]. All the models are able to correct the distortion at sub-pixel precision better than 0.1 pixel and the polynomial and the rational models seem to perform better when the distortion is large. The polynomial model raises less numerical risk since it can be solved linearly. The only requirement is to dispose of many more control points than the number of parameters. This might be an objection when using over-simple calibration patterns. In our experiments the number of control points (about 4000) was far higher, about 60 times the number of polynomial coefficients ((7 + 1) • (7 + 2) = 72 for a 7th-order polynomial). Clearly, this also entails two methodological changes: first that the distortion should be corrected independently and previously to projective calibration. Second, that calibration patterns should contain a higher than usual number of control points (more than 500). With the current image resolution in most cameras, this is no longer a restriction. With these provisions, the polynomial model is clearly the best model in terms of high precision and of ease of evaluation.

  24mm f/4G ED-IF AF-S DX ptlens Nikkor 16-35mm f/4G ED-AFS VR ptlens Nikkor 16-85mm f/3.5-5.6G AF-S ED VR DX ptlens Sony Sony AF DT 16-105mm F3.5-5.6 poly5 Sony DT 18-55mm F3.5-5.6 SAM SAL 1855 poly3 Minolta/Sony AF DT 18-70mm F3.5-5.6 50mm f/2.8 XR Di II LD ptlens Tamron 18-200mm f/3.5-6.3 XR Di II LD ptlens Pentax SMC Pentax DA 12-24mm F/4 ED AL IF poly3 SMC Pentax DA 18-55mm f/3.5-5.6 ptlens SMC Pentax DA 50-200mm f/4-5.6 DA ED ptlens

M

  with Da and Dc the energy computed on the other M pairs of points, which were not used in the energy minimization. We used points on regular grids in [-1, 1] 2 with M = 20 × 20 = 400.

Fig. 3 :

 3 Fig. 3: The geometric constraint to correct the distortion synthesized by the radial+tangential model. (a) The correction model is the radial+tangential model. The distortion center is fixed and known at (x c , y c ). (b) The correction model is the radially symmetric. The distortion center (x c , y c ) is also estimated.

Fig. 2 :

 2 Fig.2: Comparison of infinitely differentiable radial and division models with few parameters, used for correction. From left to right: one, two, and three parameters in the models. The plain curve represents the average residual error and the dashed curve the maximum error.

Fig. 4 :

 4 Fig. 4: (a) digital texture pattern. (b) and (c) two similar photographs of the flat pattern.

Fig. 5 :

 5 Fig. 5: The experiment with cameras Canon EOS 30D with EF-S 18-55mm lens. The focal length is 18mm. Top row: distorted image. Second row: corrected images. Third row: distortion field coded as the module of the displacement vector pointing from the digital pattern to the photo. The radial model and the radial+tangential model are with H and S in[START_REF] Rudakova | Camera matrix calibration using circular control points and separate correction of the geometric distortion field[END_REF]. In fact, the results of all the tested parametric models are visually identical, according to the accuracy shown in TablesIV and V.

TABLE I :

 I Some lenses in LensFun used for synthetic tests. Note that our experiments are run on a complete list of lenses available in LensFun. Refer to LensFun website for more information. The distortion models referred here are written in Table

TABLE III :

 III The synthesized distortion by one radial+tangential (R+T) model, corrected by different models. The average residual error is represented in pixels in an image of 1761 × 1174 pixels.

		-2									-2									-2
	average/max residual error (in log_10)	-16 -14 -12 -10 -8 -6 -4					radial 1 parameter division 1 parameter	average/max residual error (in log_10)	-16 -14 -12 -10 -8 -6 -4					radial 2 parameters division 2 parameters	average/max residual error (in log_10)	-16 -14 -12 -10 -8 -6 -4					radial 3 parameters division 3 parameters
		0	500	1000	1500	2000	2500	3000	3500	4000	0	500	1000	1500	2000	2500	3000	3500	4000	0	500	1000	1500	2000	2500	3000	3500	4000
						lens index									lens index									lens index

IV. REAL DISTORTION FITTING EXPERIMENTS

We now present real tests to compare the performance of the models for real distortion correction. Our test was based on a non-parametric camera distortion estimation method

[START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF]

TABLE IV

 IV

	: The average fitting error for the considered models
	for the experiment with cameras Canon EOS 30D and Canon
	EOS 40D. R+T: radial + tangential model, R: radial model,
	D: division model, F: FOV model, P: polynomial model and
	Ra: rational model.

TABLE V :

 V The average straightness error for different models for the experiment with cameras Canon EOS 30D and Canon EOS 40D. R+T: radial + tangential model, R: radial model, D: division model, F: FOV model, P: polynomial model and Ra: rational model.

To be exact, it depends on the concavity or convexity (or their absence) of the polynomial on the right hand side of (3) over the extent of the image.

This formula is slightly different from the original one proposed in[START_REF] Devernay | Straight lines have to be straight[END_REF] because we work in the normalized image domain and we consider ω as half of the field of view instead of the full field of view.

Originally the authors proposed to complete the FOV model with evenorder coefficients

k 4 , k 6 , • • • .[START_REF] Byrod | Fast and robust numerical solutions to minimal problems for cameras with radial distortion[END_REF] See http://lensfun.sourceforge.net/calibration/ for the lens profile calibration

procedure.[START_REF] Clarke | The principal point and CCD cameras[END_REF] Actually, the larger image dimension is mapped to [-1.0, 1.0] and the other dimension scaled so as to preserve the aspect ratio.

ACKNOWLEDGMENTS

This work was supported by Agence Nationale de la Recherche ANR-09-CORD-003 (Callisto project), Direction Générale de l'Armement ANR-12-ASTR-0035 (STEREO project), Centre National d'Etudes Spatiales (MISS project), Office of Naval research under grant N00014-97-1-0839, and European Research Council, advanced grant "Twelve labours".