R Grompone Von Gioi

P Monasse

J.-M Morel

Z Tang

LENS DISTORTION CORRECTION WITH A CALIBRATION HARP

Keywords: Lens distortion, error compensation, plumb lines, polynomial model

Plumb line lens distortion correction methods permit to avoid numerical compensation between the camera internal and external parameters in global calibration method. Once the distortion has been corrected by a plumb line method, the camera is ensured to transform, up to the distortion precision, 3D straight lines into 2D straight lines, and therefore becomes a pinhole camera. This paper introduces a plumb line method for correcting and evaluating camera lens distortion with high precision. The evaluation criterion is defined as the average standard deviation from straightness of a set of approximately equally spaced straight strings photographed uniformly in all directions by the camera, so that their image crosses the whole camera field. The method uses an easily built "calibration harp," namely a frame on which good quality strings have been tightly stretched to ensure a very high physical straightness. Real experiments confirm that our method produces high precision corrections (less than 0.05 pixel), approximating the distortion with a large number of degrees of freedom given by a polynomial model of order eleven.

INTRODUCTION

This paper presents a method to correct camera lens distortion with high precision. By high precision, we mean deviations from straightness of about 0.1 pixel for a straight line crossing the whole camera field. Such a precision is hardly appreciable for a human observer. However, there is no limit to the desired precision when the camera is used for 3D reconstruction or photogrammetry tasks. Traditionally, lens distortion and the other camera parameters are estimated simultaneously as camera internal and external parameters [START_REF] Slama | Manual of Photogrammetry[END_REF][START_REF] Roger | A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-theshelf 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]. In these global calibration methods all parameters are estimated by minimizing the error between the camera and its numerical model on feature points identified in several views, all in a single non-linear optimization. The result will be precise if (and only if) the model captures the correct physical properties of cameras and if the minimization algorithm finds a global minimum. Unfortunately, global camera calibration suffers a common drawback: errors in the external and internal camera parameter can be compensated by opposite errors in the distortion model. Thus the residual error can be apparently small, while the distortion model is not precisely estimated [START_REF] Weng | Camera calibration with distortion models and accuracy evaluation[END_REF][START_REF] Li | A non-iterative method for correcting lens distortion from nine point correspondences[END_REF]. For example, the Lavest et al. method [START_REF] Lavest | Do we really need accurate calibration pattern to achieve a reliable camera calibration?[END_REF] measures the non-flatness of a pattern and yields a remarkably small re-projection error of about 0.02 pixels, while the straightness of corrected lines has a 0.2 pixel RMSE. Fortunately the error compensation in global calibration can be avoided by proceeding to distortion correction before camera calibration. Recent distortion correction methods use the correspondences between two or several images, without knowledge of any camera information. The main tool they use is slackened epipolar constraints, which incorporate lens distortion into the epipolar geometry. Several iterative [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] or non-iterative methods [START_REF] Pajdla | Estimation of omnidirectional camera model from epipolar geometry[END_REF][START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF][START_REF] Barreto | Fundamental matrix for cameras with radial distortion[END_REF][START_REF] Li | A non-iterative method for correcting lens distortion from nine point correspondences[END_REF][START_REF] Claus | A rational function lens distortion model for general cameras[END_REF] are used to estimate the distortion and to correct it. These methods are used with a low order parametric distortion model and therefore cannot achieve high precision.

Non-parametric methods which establish a direct diffeomorphism between a flat pattern and a frontal photograph of it [START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF][START_REF] Stevenson | Nonparametric correction of distortion[END_REF] should be ideal for high precision distortion correction. Indeed, they do not depend on the a priori choice of a distortion model with a fixed number of parameters. Yet to achieve a high precision, they depend on the design of a very flat non deformable plate with highly accurate patterns printed on it. 1 This replaces a technological challenge by another, which is not simpler. Plumb-line methods [START_REF] Brown | Close-range camera calibration[END_REF] should therefore be an alternative because, as we shall see, it is easier to create very straight lines. For plumb-line methods, a distortion model is still necessary to precisely remove the distortion, and most existing models can be used. Nevertheless, some of them are too complicated [START_REF] Brown | Close-range camera calibration[END_REF], while some are not general enough to capture the distortion [START_REF] Devernay | Straight lines have to be straight[END_REF]. For most distortion models, the distortion center is a sensitive parameter when a realistic distortion is treated. The barely polynomial approximation proposed in [START_REF] Grompone Von Gioi | Self-consistency and universality of camera lens distortion models[END_REF] is therefore a good choice, being a translation invariant and linear approximation of any vector field. This model-free formulation can approximate complex radial and non-radial distortions as well, provided the polynomial degree is high enough. According to the criteria of self-consistency and universality2 developed in [START_REF] Grompone Von Gioi | Self-consistency and universality of camera lens distortion models[END_REF] to compare many camera distortion models, the polynomial models are the most flexible and accurate.

The proposed method is introduced in section 2, followed by real experiments in section 3, along with a comparison to a non-parametric method. Section 4 concludes the paper.

THE HARP CALIBRATION METHOD

In one sentence, the proposed method combines the advantage of plumb-line methods with the universality of the polynomial approximation. The plumb-line method consists in correcting the distorted points which are supposed to be on a straight line, by minimizing the average squared distance from the corrected points to their corresponding regression lines.

The polynomial model has the form

x u = b 0 x p d + b 1 x p-1 d y d + b 2 x p-2 d y 2 d + • • • + b p y p d +b p+1 x p-1 d + b p+2 x p-2 d y d + • • • + b 2p y p-1 d + • • • + b (p+1)(p+2) 2 -3 x d + b (p+1)(p+2) 2 -2 y d +b (p+1)(p+2) 2 -1 y u = c 0 x q d + c 1 x q-1 d y d + c 2 x q-2 d y 2 d + • • • + c q y q d +c q+1 x q-1 d + c q+2 x q-2 d y d + • • • + c 2q y q-1 d + • • • + c (q+1)(q+2) 2 -3 x d + c (q+1)(q+2) 2 -2 y d +c (q+1)(q+2) 2 -1 (1)
with (x u , y u ) undistorted point and (x d , y d ) distorted point. The polynomial approximation being translation invariant, the origin is arbitrarily fixed at the image center. The order for the x and y components is respectively p and q. The number of parameters for x and y is respectively (p+1)(p+2) 2 and (q+1)(q+2)

2

. The model is called bicubic model when p = q = 3.

In the following, we show how to integrate the polynomial model into the plumb-line method. Given a set of corrected points (x ui , y ui ) i=1,••• ,N which are supposed to be on a line, the first step is to compute the linear regression line 2 , and γ = A x sin θ + A y cos θ. The sum of squared distances from the points to this regression line is

αx ui + βy ui -γ = 0 (2) with α = sin θ, β = cos θ, tan 2θ = - 2(Axy-AxAy) Vxx-Vyy , A x = 1 N N i=1 x ui , A y = 1 N N i=1 y ui , A xy = 1 N N i=1 x ui y ui , V xx = 1 N N i=1 (x ui -A x ) 2 , V yy = 1 N N i=1 (y ui -A y )
N i=1 (αx ui + βy ui -γ) 2 .
By considering G groups of lines, the total sum of squared distance is

S = G g=1 Lg l=1 N gl i=1 (α g x u gli + β g y u gli -γ gl ) 2 (3)
with L g the number of lines in group g, and N gl the number of points of line l in group g. Given the total number of points

N = G g=1
Lg l=1 N gl , the root mean squared distance (RMS error) is defined by

d = S N . (4)
For a sake of succinctness, the following discussion will assume a bicubic model with p = q = 3. Combining Eq. ( 1) and Eq. ( 3), the energy S becomes

S = G g=1 Lg l=1 N gl i=1 α g b 0 x 3 d gli + • • • + b 9 ( 5
)
+β g c 0 x 3 d gli + • • • + c 9 -γ gl 2 .
It is a non-linear problem to minimize the energy S with respect to the parameters b 0 ,

• • • , b 9 , c 0 , • • • , c 9 .
But the problem becomes linear by assuming that the orientation parameters α g , β g are known. By differentiating S with respect to each parameter, we obtain a linear system In practice, the orientation of lines is unknown and the minimization of the energy in Eq. ( 5) is a non-linear problem. The minimization is performed by first doing an iterative Levenberg-Marquardt (LM) algorithm which estimates the parameters of polynomials of increasing order. The algorithm starts estimating the parameters of a 3-order polynomial; the result is used to initialize the 4-order polynomial, and the process continues until 11-order. After this first step, the linear estimation in Eq. ( 6) is performed iteratively to refine the precision. The line orientations are first initialized by the orientation of the regression lines obtained by the LM method, and then with the values of the previous linear step. The iteration is repeated until the results stabilize or the required precision is reached.

Ax = 0 (6) with x = (b 0 , • • • , b 9 , c 0 , • • • , c 9 , γ 11 , • • • , γ GL G )

EXPERIMENTS

In this section, we describe real experiments with a strong distortion comparing the proposed method with a nonparametric flat pattern-based method [START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF].

The experiments were made with a Canon EOS 30D reflex camera with EFS 18 -55mm lens. The minimal focal length (18mm) was chosen to produce a fairly large distortion. The RAW images were demosaicked by summing up the four pixels of each 2 × 2 Bayer cell, obtaining a half-size image. The calibration harp was built by tightly stretching good quality strings on a frame to ensure a very high physical straightness (our calibration harp was built with sewings string). The high distortion is visible near the border of the image (see the images in Fig. 1 for example). On each corrected line, sub-pixel precise edge points were obtained by Devernay's algorithm [START_REF] Devernay | A non-maxima suppression method for edge detection with sub-pixel accuracy[END_REF] and then groupped when belonging to the same line segments detected by the LSD algorithm [START_REF] Grompone Von Gioi | LSD: A fast line segment detector with a false detection control[END_REF]. 18 photographs of the calibration harp with different orientations (some of them in Fig. 1) were used in the calibration by the 11-degree polynomial model. Independent distorted images (for example, see the top image in Fig. 2) are used for verification. The correction result is recapitulated in Fig. 2. The non-parametric pattern-based method in [START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF] was also tried (see result in Fig. 3), which estimates the distortion as the diffeomorphism (up to a homography) mapping the original digital pattern to a photograph of it by triangulating and interpolating dense correspondences. Note that in Fig. 2 and3, y-axis has the same range, from -0.3 pixels to 0.3 pixels.

For the non-parametric method, a global tendency in the straightness error of the corrected lines can be observed (see Fig. 3). This was in fact due to the unavoidable drawback of this method: there is never a guarantee that the pattern is completely flat. The non-flatness of the pattern introduces a bias in the estimated distortion field, which causes the observable global residual in the curves plotted in Fig. 3. To eliminate this error source, the solution is either to construct a very flat pattern, or to recover the 3D shape of a non-flat pattern. Both are not really feasible. Yet, to appropriately use a plumb-line method, we need a pattern containing very straight lines, and this is far easier in practice. As shown in Fig. 2, the distortion correction is so accurate that no visible global tendency is visible in the corrected curves. The root mean square (RMS) distance of each line is also significantly smaller than for the non-parametric method (Table 1). It is particularly striking in Fig. 2 that the superimposed curves of the left and right side of each string are fairly uncorrelated, meaning that no deterministic distortion is left. The erratic oscillation of very small amplitude can be attributed to any cause, from the lack of the uniformity of the harp background causing a shift in the edge detection, to aliasing in the image itself. But it cannot be due to a residual mismatch of the polynomial model itself, because otherwise the curves on both sides of each string would be parallel. This confirms a posteriori the reliability of the polynomial model. RMSE (in pixels) proposed method non-parametric [ 1: RMS Error from edge points of corrected lines to their corresponding regression line. The proposed method is compared to non-parametric pattern-based method [START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF]. Each cell in the table contains two values, one for each side of the string. The lines are numerated in the top image in Fig. 2.

explains why we can call the method model-free. The only assumption on the lens distortion is its smoothness, implying that a polynomial with high enough order approximates it. In our experiments, the approximation error stabilizes for polynomials of degree 7 to 11. It might be objected that the high number of parameters in the polynomial interpolation (156 for an 11-order polynomial) could cause some bias in the result. Yet, the number of control points is far higher: There were about 10 strings for each orientation, some 30 control points on each string side, and some 18 orientations. Thus the number of control points is about 10000 and therefore 60 times more than the number of polynomial coefficients. A visual examination of the two sides of the strings confirms that no artificial simultaneous bias has been introduced by the polynomial distortion correction. This observation seems to indicate that most of the 0.05 pixels remaining oscillation is due either to image processing factors, or to background inhomogeneity, to aliasing in the edge detector, or to string diameter variations. It is not easy to decide which factor is dominant. Future work will concentrate on the precision of external camera calibration, and eventually on the 3D precision after having removed the lens distortion by the present method.

  T and the coefficient matrix A. To obtain a unique non-trivial solution, we always set b 7 = c 8 = 1, b 9 = c 9 = 0, which in fact fixes a scale and a translation to the solution. The minimized S can be changed by the introduced scale. But this change is consistent if the distortion center and b 9 , c 9 , b 7 , c 8 are fixed.

  By combining the advantages of a model-free polynomial approximation and of a real plumb line pattern, the proposed lens distortion correction is significantly more accurate than parametric methods on flat patterns. The "calibration harp" construction only requires the acquisition of a string with decent quality. It is far simpler than realizing a flat plate with highly accurate patterns engraved on it. (The calibration of such patterns is not easier than lens calibration itself!) The high number of degrees of freedom in the unstructured model line No.

Fig. 1 :

 1 Fig. 1: Six of the 18 photos of "calibration harp" with different orientation.

Fig. 2 :

 2 Fig. 2: Correction performance of the proposed method. Top row, the independent distorted image. The lines are numerated in red. From the second row to the last row, from left to right: the distance from the edge point of corrected lines to their corresponding regression line. Note that each figure contains two curves because there are two lines for one string. The range of y-axis is from -0.3 pixels to 0.3 pixels.

Fig. 3 :

 3 Fig.3: Correction performance of non-parametric patternbased method[START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF]. See the caption of Fig.2for details.

A 10 micron flatness could be needed to achieve a precision of 0.01 pixels.

Self-consistency is evaluated by the residual error when distortion generated with a certain model is corrected (using the model in reverse way) by the best parameters for the same model. Analogously, universality is measured by the residual error when a model is used to correct distortions generated by a family of other models. A model is self-consistent and universal if it can approximate any other model and the inverse of any other model, including itself, with the desired precision. Polynomials of order 11 are 0.01 pixels self-consistent and universal.


Zhongwei Tang
email: tang@cmla.ens-cachan.fr

Rafael Grompone Von Gioi

Pascal Monasse

Jean-Michel Morel

Self-consistency and universality of camera distortion models

This paper introduces the concepts of "selfconsistency" and "universality" to evaluate high precision camera distortion models. Self-consistency is evaluated by the residual error when the distortion generated with a certain model is corrected by the best parameters for the same model (used in reverse way, which is common practice). Analogously, universality is measured by the residual error when a model is used to correct distortions generated by a family of other models. Five classic camera distortion models are reviewed and compared for their degree of self-consistency and universality. Among the evaluated models, it is concluded that only the polynomial and the rational models are universal up to precisions of 1/100 pixel. However, the polynomial model, being linear, is much simpler and faster to estimate. Unusually high polynomial degrees are required to reach this strong precision. Nevertheless, extensive numerical experiments show that such distortion polynomials are easily estimated and produce a precise distortion correction without over-fitting. Our conclusions are validated by three independent experimental setups: The models are compared first in synthetic experiments by their approximation power; second by fitting a real camera distortion estimated by a non parametric algorithm; and finally by the absolute correction measurement provided by photographs of tightly stretched strings, warranting a high straightness.

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, mainly because of lens geometric distortion [START_REF] Brown | Close-range camera calibration[END_REF], and possibly from the CCD shape distortion itself. 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 computing power, high-precision 3D reconstructions become feasible, demanding 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 to three orders of magnitude. This increased accuracy requires a new methodology for evaluating distortion models. In a nutshell, our conclusion is that a polynomial model of higher degree than usual, ranging from 8 to 15, is necessary for reaching a pixel precision ranging from 1/100 to 1/1000. The polynomial model permits to approximate at this resolution any other model, and the inverse of any other model, including itself. When these properties are reached, the model is called universal and self-consistent. Among the other four models which will be compared (radial, division, FOV, and rational), only the rational model has the eligible self-consistency and universality, but to a far higher computational cost.

With the exception of a few non-parametric methods [START_REF] Farid | Blind removal of lens distortion[END_REF], [START_REF] Stevenson | Nonparametric correction of distortion[END_REF], [START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF], an appropriate distortion model is indispensable to model the deviation of a real camera from an ideal pinhole camera. The main distortion models are the radial model [START_REF] Brown | Close-range camera calibration[END_REF], the division model [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF], the FOV model [START_REF] Devernay | Straight lines have to be straight[END_REF], the bicubic model [START_REF] Kilpelä | Compensation of systematic errors of image and model coordinates[END_REF], the rational model [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]. This diversity is only marginally linked to the kind of camera. Thus, a synthetic quantitative and qualitative comparison is required. Do these models reflect camera distortion in its physical aspect? It could be argued that a correct model should originate from physical measurements on systems of lenses. Surprisingly enough, there is little physical background for the distortion models in the literature. In [START_REF] Weng | Camera calibration with distortion models and accuracy evaluation[END_REF], lens distortion is decomposed into three effects: radial distortion, decentering distortion and thin prism distortion. But, still, it is only marginally based on a physical background. In fact, the final distortion includes effects caused by a complex lens system, by the camera geometry, and by the (not perfectly planar) shape of the captor. One is therefore led to figure out a flexible model with enough parameters to approximate any plausible distortion. In absence of a physical model, the model classification approach adopted here will be to look for models which actually cope with any other proposed distortion model, at a given precision.

The second question is the relationship between the distortion and the correction models, which should be inverse of each other. A distortion model is used to simulate the distortion of ideal images, while a correction model is used to correct distorted images. Indeed, most of the widely used models are not invertible, thus the correction model and the distortion model must be different. In the literature, however, it seems that the roles of distorted point and undistorted point are interchangeable, which again confirms the lack of physical meaning for these models. For example, direct distortion models are used in global camera calibration [START_REF] Roger | A versatile camera calibration technique for highaccuracy 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] Gomez | Algebraic lens distortion model estimation[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 correction models are used without any fuss to approximate the inverse distortion.

Assume we simulate a camera distortion with a certain model and a certain set of parameters. Except for some models, the distortion will not be corrected by using the same model with other parameters, because the model itself is usually not invertible. We propose to measure the error incurring when a distortion is inverted by using the same model as the one used for simulating the distortion. This inversion error, when the best correcting parameters have been estimated and applied, will be a measurement of the model self-consistency. In other words, self-consistency relates to how well a model is able to correct distortion generated by a model of its own family. Of course the best models should be universal, therefore able to correct distortions generated by other models. We therefore propose to measure a model universality as the residual error when this model is used to correct distortion generated by a whole set of different models. A universal model is a model for which this error is very small no matter what other (reasonable) distortion model has been applied. A selfconsistent and universal model implies that the model can also approximate the distortion generated by other models. Our goal is to identify the least complex universal and selfconsistent models. Of course this question only makes sense within fixed accuracy bounds. None of the classic models is actually exactly (algebraically) invertible. Thus, the introduced concepts, universality and self-consistency, must be thought of with the to a given precision caveat. As a matter of fact, for off-the-shelf cameras, most distortion models are roughly equivalent at a 1 pixel precision. The question is different when we aim at sub-pixel precisions. These precisions, up to 1/100 pixel, are highly desirable when using cameras for stereovision or photogrammetric tasks.

The other caveat is that, although distortion models reflect a model of the optical lens, the real corrected distortion must actually involve the whole system lens + 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 realistic synthetic distortion data permitting to quantify the ideal attainable precision. Then, the same models will be compared on their capacity to fit a real camera distortion (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 [26]. 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 self-consistency and universality are evaluated in Section III by synthetic experiments. Section IV and V describe the experiments done with real cameras. Section VI is a conclusion.

II. DISTORTION AND CORRECTION MODELS

Let us denote by (x u , y u ) the coordinates of an undistorted point as would be observed in an ideal pinhole camera. Due to the lens geometric distortion, this point will be distorted to coordinates (x d , y d ). We will model distortion 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 transformation in the opposite direction:

x u = g x (x d , y d ), y u = g y (x d , y d ). (2)
A particularly interesting case is when the functions f or g show radial symmetry relative to a fixed distortion center (x c , y c ). In that case we obtain a compact formulation using normalized 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, namely the radial model [START_REF] Brown | Close-range camera calibration[END_REF], the division model [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF], the FOV model [START_REF] Devernay | Straight lines have to be straight[END_REF], the polynomial model [START_REF] Kilpelä | Compensation of systematic errors of image and model coordinates[END_REF], and the rational function model [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]. 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 distorted 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 representing a scaling does not introduce distortion. The scaled image is distorted by k Note that the distortion center (x c , y c ) is also a parameter of radial models. The division model is nothing but the scalar inverse of the radial model,

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 ω:

r 2 = r 1 tan(r 1 ω) 2r 1 tan( ω 2 ) . ( 5
)
Here the coefficient ω is of order 1, but more coefficients can be added to the FOV model to make it more complete.

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 . (6)
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 . (7)

III. SELF-CONSISTENCY AND UNIVERSALITY

In the literature it is not always clear whether the above models are correction models or distortion models. This raises the question of self-consistency, while the universality question is raised by the plurality of models. Being theoretical properties of model families, both properties can be genuinely evaluated by synthetic experiments. Self-consistency and universality will be tested by generating a distortion with any of the above models, and then evaluating the residual error after applying the best distortion correction with each model. Of course, we will use for each model sets of parameters that generate a realistic distortion, as illustrated in Table I. A distortion is generated using a first model in direct way, using Eq. ( 1). This distortion is then corrected by identifying the best parameters of a second model when used as in Eq. ( 2).

In our synthetic tests, both (x u , y u ) and (x d , y d ) are known. The question is 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 g x and g y by minimizing the difference between the ideal correction and the practical correction. The energy to be minimized can be written as

C = g x (x d , y d ) -x u 2 + g y (x d , y d ) -y u 2 dx d dy d . (8)
In practice, the simulation is performed on M samples (x ui , y ui ), i = 1, . . . , M regularly distributed on an image. The corresponding distorted samples (x di , y di ), i = 1, . . . , M are obtained by Eq [START_REF] Alvarez | An algebraic approach to lens distortion by line rectification[END_REF]. The discrete energy to be minimized is

D = M i=1 g x (x di , y di ) -x ui 2 + g y (x di , y di ) -y ui 2 . ( 9
)
Given the two coordinates for each of the M points, we have 2M equations relating the ideal and the distorted points. In the case of the polynomial model these equations have the form given in Eq. ( 6), and the problem can be formulated as a linear system Ak = b, [START_REF] Grompone Von Gioi | LSD: A fast Line Segment Detector with a false detection control[END_REF] where A is a matrix containing the different powers for all the distorted points, k is formed with the model's coefficients, and b contains the undistorted coordinates. The set of model coefficients with least error D is obtained by minimizing the norm Akb 2 , which results in

k = A T A -1 A T b. (11)
In practice, the matrix A is ill-conditioned and can make the solution unstable. Some normalization technique should be applied before.

The same linear method can be applied to the radial model when the distortion center is known because the transformation Eq. ( 3) is also linear. 1 For all the other models, a non-linear method must be used, even if (x c , y c ) is known. The minimization is performed by first doing an incremental Levenberg-Marquardt (LM) algorithm which estimates the parameters in increasing order. The algorithm starts estimating the parameters of a low 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 is computed explicitly to make the algorithm efficient. Even though this strategy is complex, it avoids some local minima and is safer than performing LM directly on the model of the aimed order.

A. Experiments with known distortion center

For our first experiment, a total of M = 5104 points were regularly distributed in an image domain of size 1761 × 1174 and the distortion center was fixed at (880.5, 587) and was assumed to be known. The corresponding distorted points were computed using the different models with the parameters in Table I. The self-consistency and universality measurements are recapitulated in Table II, expressed as the average error D = D M after estimating the parameters which minimize the energy in Eq [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF].

B. Experiments with unknown distortion center

In practice, the distortion center (x c , y c ) is unknown and it should also be considered as a parameter in the minimization formulation. Our second experiment was done in this way. The minimization problem becomes now non-linear for most models. This is true for the radial model, the division model and the FOV model. In contrast, the polynomial and the rational function models are invariant to a translation of the distortion center. The point (x c , y c ) can be fixed arbitrarily, and in the polynomial case the minimization problem is still linear. This is a decisive advantage with respect to the other models. The self-consistency and universality results are recapitulated in Table III with the same parameters for generating distortion in Table I. For the distortion generation, the distortion center was fixed again at the center (880.5, 587) of the image, while for the correction the initial distortion center was realistically taken (50, 50) pixels away from the true position. For the radial model and the division model, the Levenberg-Marquardt algorithm could still find the true distortion center, and the minimized error was the same as when the distortion center was known. Nevertheless, for the FOV model, a wrong initialization of the distortion center degraded the correction performance. For the polynomial model, the solution can be found linearly by fixing an arbitrary distortion center.

Distortion Generation Model R 2 • R 4 • D 2 • D 4 • F 3 • P 3 • P 4 • Ra 2 • Correction Model R 2 • 0.

C. Comparison

We aim at an average precision below 10 -2 pixel. The tables show that the models are self-consistent if the order of correction is high enough. The radial model and the division model are consistent with each other, whether the distortion center is known or not. The FOV model is a little less consistent with the radial symmetric models, including itself, when the distortion center is known. With an unknown distortion center, the FOV correction performance decays. The polynomial model instead seems to be able to correct any type of distortion, but a higher order is often necessary to correct the radial, division or FOV distortions. This higher order is not a problem, because of the computational efficiency of the linear method. The rational model has a comparable performance at the price of a much higher computational cost. In conclusion, at precision 1/100 pixels the polynomial and the rational models are the only ones to be jointly selfconsistent, universal and linear among the compared models. But the polynomial model has a much lower complexity than the rational model to attain the same precision. A degree 8 polynomial model is enough to attain precisions strictly below 1/10. The 1/100 pixel precision is robustly attainable with degree 15 polynomials.

D. Realistic distortion

A real distortion can be far more complex than what the above simple models can generate. A more realistic distortion contains a radial symmetric term, a term for decentering distortion and a term for thin prism distortion [START_REF] Weng | Camera calibration with distortion models and accuracy evaluation[END_REF], 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 radial symmetric distortion and tangential distortion. In Table IV, the self-consistency and universality of the models were again tested with known distortion center after adding the additional distortion with p 1 = 4.0 × 10 -6 , p 2 = -2.0 × 10 -6 , p 3 = 0, s 1 = 3.0 × 10 -6 , s 2 = 1.0 × 10 -6 . By adding a non-radial component in the distortion, the radial model, the division model and the FOV model do not reach anymore the 10 -2 pixel precision. Both polynomial model and rational function model give a precision equal or better than 10 -2 pixel with a high model order. But it is always linear to solve the parameters of polynomial model, while it requires a complex incremental LM minimization to solve rational function model.

xd = xu k 0 + k 1 r u + k 2 r 2 u + • • • + p 1 r 2 u + 2x 2 u + 2p 2 xu ȳu 1 + p 3 r 2 u + s 1 r 2 u , ȳd = ȳu k 0 + k 1 r u + k 2 r 2 u + • • • + p 2 r 2 u + 2ȳ 2 u + 2p 1 xu ȳu 1 + p 3 r 2 u + s 2 r 2 u , Distortion Generation Model R 2 • R 4 • D 2 • D 4 • F 3 • P 3 • P 4 • Ra 2 • Correction Model R 2 • 0.

IV. REAL DISTORTION FITTING EXPERIMENTS

After its validation on synthetic examples, we present here real tests to verify that the proposed high order polynomial model works for real distortion correction. This test is based on the non-parametric camera distortion estimation method in [START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF] 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 mirror was used in the experiments). Two photos of the pattern were taken by a Canon EOS 30D SLR camera with EFS 18-55mm lens. The minimal focal length (18mm) was chosen (with fixed focus) to produce a fairly large distortion. The distortion was estimated (up to a homography) as the diffeomorphism mapping the original digital pattern to a photograph of it. Without going into details, the algorithm is summarized in the following:

1) Take two slightly different photographs of a textured planar pattern with a camera whose settings are frozen; 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 both photographs to obtain matching pairs of points; 3) eliminate outliers by a loop validation step; 4) triangulate and interpolate the remaining matches to get a dense reverse distortion field; 5) 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; 6) by applying the reverse distortion field to all images produced by the same camera, the camera is converted into a virtual pinhole camera. The matching pairs delivered by step 5 (about 8000 in our experiments) in the above algorithm are "outliers"-free and precise thanks to the loop validation and filtering by local homography. So we can directly try all models to fit these "outliers"-free 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 classically 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   . (12)
Since the polynomial model and the rational function model can approximate well H and S, we can use these models to approximate the distortion field without explicitly estimating the homography. Nevertheless, for the radial symmetric models, it is necessary to take into account H and S when approximating the distortion. Indeed H and S are generally not radial symmetric. Thus, we have 9 additional parameters to estimate, besides the parameters of radial symmetric distortion model. 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 and the homography.

Half of the matching pairs were used to estimate the parameters for the different models, and the other half to evaluate the average fitting error. The results are recapitulated in Table VI. They show that by combining H and S to model the inclination between the camera and the pattern, all of the radial symmetric models give almost the same fitting error, which becomes stable (0.15 pixel) when the order attains 4. Table V shows H and S estimated by three radial symmetric models of order 12. The similar estimation of H and S implies that the minimization process is stable.

Distortion Generation Model R 2 • R 4 • D 2 • D 4 • F 3 • P 3 • P 4 • Ra 2 • Correction Model R 2 • 7
The polynomial model and the rational function model give a stable fitting error (0.04 pixel), which is about 4 times smaller than the radial 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 fitting error of the rational function model becomes stable when its order attains 4 (45 parameters to be estimated), while the polynomial model gives a stable fitting error when its order attains 7 (72 parameters to be estimated). Even though in this experiment the rational function model converges faster than the polynomial model, solving the rational function model requires a complex incremental LM minimization, which is a time-consuming process requiring a good initialization and does not always ensure the convergence to the global minimum, 2 while the polynomial model can always be solved linearly.

The fact that the polynomial model and the rational function model give fitting errors about 4 times smaller than the radial symmetric models implies that the real camera distortion is not strictly radial symmetric. Otherwise, the fitting error of the radial symmetric models should be at least as small as 0.04 pixel. The residual fitting error (0.04 pixel) given by the polynomial model and the rational function model can be attributed to the noise of matching points.

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 need a more objective evaluation to check the quality of the correction models. To this purpose, a physical frame with tightly stretched cylindrical strings was built [26]. 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 2 The linear solution to the rational function model proposed in [START_REF] Claus | A plumbline constraint for the rational function lens distortion model[END_REF] images of strings taken by the same camera with the same fixed lens configuration (see Fig. 1). The average distance from the edge points (computed by the method of [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 was computed. See [26] for more details. 12). Third row: the correction result by the polynomial model of order 12. From left to right column: the distortion field coded as the module of the displacement vector pointing from the digital pattern to the photo; the level line on the distortion field with quantification step of 10; corrected image. The results of the radial model and the polynomial model are visually identical. In fact, the results of all the tested parametric models are visually identical, according to the accuracy shown in Table VI and VII. The non-parametric method gives a slightly different result at the border due to the fact that the triangulation is imprecise at the border. models are not capable to capture.

VI. CONCLUSION

We introduced the self-consistency and universality criteria for camera distortion models. Using these tools, five classic distortion models were evaluated and compared. The polynomial and rational model were shown to be both selfconsistent and universal, to the cost of a high degree. This high degree raises no computational issue for the polynomial model. Indeed, we have seen that after a correct conditioning it can always be solved linearly (in contrast to the rational model). Furthermore, the polynomial model is translation invariant, which makes it insensitive to a translation of the distortion center. This model is not adapted to global camera calibration methods where the internal and 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. It might be objected that the high number of parameters in the polynomial interpolation (156 for an 11-order polynomial) could cause over-fitting bias in the results. 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 30 times the number of polynomial coefficients. Our experiments show that the residual errors stabilize for orders between 6 to 12, confirming that no over-fitting occurred. Our experiments also show that [START_REF] Grompone Von Gioi | LSD: a Line Segment Detector[END_REF].

there is a small component of non radial symmetric distortion introduced by the camera, which cannot be modeled by radial symmetric models. High order polynomials are really needed to capture the non radial symmetric distortion if we wish to obtain high precisions. 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 more a restriction.

model parameters radial 2 •

 2 k 0 = 1.0, k 1 = 0.25 × 10 -4 , k 2 = -0.5 × 10 -7 1084 → 1050 radial 4 • k 0 = 1.0, k 1 = 0.25 × 10 -4 , k 2 = -0.5 × 10 -7 , k 3 = 1.0 × 10 -10 , k 4 = -1.5 × 10 -14 991.6 → 1050 division 2 • d 0 = 1.0, d 1 = -0.25 × 10 -4 , d 2 = 0.5 × 10 -7 1083 → 1050 division 4 • d 0 = 1.0, d 1 = -0.25 × 10 -4 , d 2 = 0.5 × 10 -7 , d 3 = -1.0 × 10 -10 , d 4 = 1.5 × 10 -14 988.7 → 1050 FOV 3 • k 0 = 1.0, ω = 1.0 × 10 -3 , k 2 = -2.0 × 10 -7 , k 3 = 4.0 × 10 -10 501.4 → 1050 polynomial 3 • a 1 = b 1 = -1.0 × 10 -8 , . . ., a 5 = b 5 = 2.0 × 10 -5 , . . ., 1050 → 1064 a 8 = 0.9, a 9 = 0.1, a 10 = 0.0, b 8 = 0.1, b 9 = 0.9, b 10 = 0.0 polynomial 4 • a 1 = b 1 = 5.0 × 10 -12 , a 6 = b 6 = -1.0 × 10 -8 , a 10 = b 10 = 2.0 × 10 -5 , . . ., 1050 → 1075 a 13 = 0.9, a 14 = 0.1, a 15 = 0.0, b 13 = 0.1, b 14 = 0.9, b 15 = 0.0 rational 2 • a 1 = 1.0 × 10 -5 , a 2 = 2.0 × 10 -5 , a 3 = 3.0 × 10 -5 , a 4 = 0.9, a 5 = 0.1, a 6 = 0.0 1031 → 1104 b 1 = 3.0 × 10 -5 , b 2 = 2.0 × 10 -5 , b 3 = 1.0 × 10 -5 , b 4 = 0.1, b 5 = 0.9, b 6 = 0.0 c 1 = 1.0 × 10 -8 , c 2 = 1.0 × 10 -8 , c 3 = 1.0 × 10 -8 , c 4 = 0.0001, c 5 = 0.0001, c 6 = 1.0

Fig. 1 .

 1 Fig. 1. Top row: distorted image. Second row: the correction result by the radial model of order 12 with H and S in Eq.[START_REF] Grompone Von Gioi | Towards high-precision lens distortion correction[END_REF]. Third row: the correction result by the polynomial model of order 12. From left to right column: the distortion field coded as the module of the displacement vector pointing from the digital pattern to the photo; the level line on the distortion field with quantification step of 10; corrected image. The results of the radial model and the polynomial model are visually identical. In fact, the results of all the tested parametric models are visually identical, according to the accuracy shown in TableVIand VII. The non-parametric method gives a slightly different result at the border due to the fact that the triangulation is imprecise at the border.

  1 , k 2 , . . .. If k 1 , k 2 , . . . are all positive, we have a pincushion distortion; if k 1 , k 2 , . . . are all negative, a barrel distortion. Mustache distortion occurs if the signs of k 1 , k 2 , . . . are not the same.

TABLE I MODELS

 I USED TO GENERATE DISTORTION, WITH THEIR REALISTIC PARAMETERS. THE VALUES ON THE LEFT AND ON THE RIGHT OF THE SIGN "→" ARE THE UNDISTORTED RADIUS AND THE DISTORTED RADIUS RESPECTIVELY. FOR THE POLYNOMIAL MODEL, THE COEFFICIENTS ARE THE SAME FOR x AND y COMPONENT, EXCEPT FOR THE ORDER 1 COEFFICIENTS. NOTE THAT THE DISTORTION CAN BE BARREL, PINCUSHION OR MUSTACHE.

TABLE II SELF

 II -CONSISTENCY AND UNIVERSALITY WITH KNOWN DISTORTION CENTER. THE AVERAGE ERROR ( D) (IN PIXELS) IS SHOWN. THE LEFT COLUMN ENTRIES SHOW THE MODEL AND THE ORDER USED FOR CORRECTION. THE TOP ENTRY ROW GIVES THE MODEL AND THE ORDER USED TO GENERATE THE DISTORTION. THE FIVE COMPARED MODEL CLASSES ARE R-RADIAL, D-DIVISION, F-FOV, P-POLYNOMIAL, AND RA-RATIONAL. THE PARAMETERS IN TABLE I WERE USED TO GENERATE THE DISTORTION. THE BOLD FONT IS USED TO HIGHLIGHT THE AVERAGE ERROR D 10 -2 .

TABLE III SELF

 III -CONSISTENCY AND UNIVERSALITY WITH UNKNOWN DISTORTION CENTER. THE INITIAL DISTORTION CENTER WAS SET (50, 50) PIXELS AWAY FROM ITS TRUE POSITION. THE AVERAGE ERROR ( D) (IN PIXELS) IS SHOWN. THE LEFT COLUMN ENTRIES GIVE THE MODEL AND THE ORDER USED FOR CORRECTION. THE TOP ROW ENTRIES GIVE THE MODEL AND THE ORDER USED TO GENERATE THE DISTORTION. THE FIVE COMPARED MODEL CLASSES ARE R-RADIAL, D-DIVISION, F-FOV, P-POLYNOMIAL, AND RA-RATIONAL. THE PARAMETERS INTABLE I WERE USED TO GENERATE THE DISTORTION. THE BOLD FONT IS USED TO HIGHLIGHT THE AVERAGE ERROR D 10 -2 .

  TABLE IV SELF-CONSISTENCY AND UNIVERSALITY WITH KNOWN DISTORTION CENTER. ADDITIONAL DECENTERING AND THIN-PRISM DISTORTION IS ADDED. COMPARE WITH TABLES II AND III. EACH ENTRY SHOWS THE AVERAGE ERROR ( D) (IN PIXELS). THE PARAMETERS IN TABLE I WERE USED TO GENERATE THE DISTORTION. THE BOLD FONT IS USED TO HIGHLIGHT THE AVERAGE ERROR D 10 -2 . THE ONLY NON-BOLD TO BOLD LINES ARE OBTAINED FOR THE POLYNOMIAL MODEL WITH DEGREE 8 TO 15.

  does not give a precise result in this case.

	order Radial Division FOV	Polynomial Rational
	3	0.29	0.29	0.28	1.48	0.19
	4	0.15	0.15	0.15	1.26	0.05
	5	0.15	0.15	0.15	0.21	0.05
	6	0.15	0.15	0.15	0.08	0.05
	7	0.15	0.15	0.15	0.04	0.04
	8	0.15	0.15	0.15	0.04	0.04
	9	0.15	0.15	0.15	0.04	0.04
	10	0.15	0.15	0.15	0.04	0.04
	11	0.15	0.15	0.15	0.04	0.04
	12	0.15	0.15	0.15	0.04	0.04
				TABLE VI
	THE FITTING ERROR (IN PIXELS) OF THE COMPARED MODELS TO THE
	MATCHINGS BETWEEN A DIGITAL TEXTURED IMAGE AND ITS
	PHOTOGRAPH (WITH HOMOGRAPHY INITIALIZED BY 8-POINT
	ALGORITHM). THE MATCHINGS OBTAINED AT STEP 5 IN THE
	SUMMARIZED ALGORITHM ARE "OUTLIERS"-FREE AND PRECISE.
	COLUMN 1 IS THE ORDER OF THE MODEL. HALF OF THE MATCHINGS ARE
	USED TO ESTIMATE THE PARAMETERS AND THE AVERAGE FITTING ERROR
	SHOWN IS COMPUTED ON THE OTHER HALF BY APPLYING THE ESTIMATED
			PARAMETERS.

  Table VII recapitulates average distance for all lines in the image. According to TableVI, the stabilized fitting error implies that none of the models has the problem of noise fitting, which guarantees the correction quality and stability. The polynomial model gives a stable performance when the order attains 6, which means that a polynomial model of order 6 was already capable of capturing the whole distortion. The rational model gives a comparable performance when its order attains 4, at the price of more complex non-linear incremental minimization which does not ensure a global minimum. The residual straightness error of corrected lines is due to the fact that the mirror on which we pasted the pattern is not completely flat. All the radial symmetric models give a straightness error larger than the polynomial model and the rational function model, due to the flatness error of the pattern and to the non radial symmetric component in the camera distortion which the radial symmetric

	model	order h 11	h 12	h 13	h 21	h 22	h 23	h 31	h 32	s
	Radial	12	1.22	0.10	-170.16	0.015	1.24	-77.76	4.8 × 10 -5	1.1 × 10 -4	1.00
	Division	12	1.22	0.10	-170.33	0.015	1.24	-77.49	4.8 × 10 -5	1.1 × 10 -4	1.00
	FOV	12	1.22	0.10	-170.15	0.015	1.24	-77.76	4.8 × 10 -5	1.1 × 10 -4	1.01

TABLE V THE

 V ESTIMATED H AND S COMBINED WITH RADIAL SYMMETRIC MODELS OF ORDER 12 TO APPROXIMATE THE DISTORTION FIELD FROM DIGITAL PATTERN TO THE PHOTO.

TABLE VII THE

 VII AVERAGED DISTANCE (IN PIXELS) FROM EDGE POINTS OF CORRECTED LINES TO THE CORRESPONDING REGRESSION LINE. THE PARAMETERS OF THE MODELS ARE ESTIMATED USING HALF OF THE MATCHINGS COMING FROM STEP 5 IN THE SUMMARIZED ALGORITHM. THE DISTORTED IMAGE IN FIG. 1 IS THEN CORRECTED BY USING ALL MODELS. THE CORRECTED LINES ARE EXTRACTED BY USING THE LINE SEGMENT DETECTOR IN [10], WHICH IS DESCRIBED IN DETAIL AND CAN BE TESTED ON LINE AT

	order	Radial Division FOV	Polynomial Rational
	3	0.27	0.27	0.27	0.74	0.11
	4	0.19	0.21	0.19	0.58	0.09
	5	0.19	0.21	0.19	0.15	0.09
	6	0.19	0.21	0.19	0.09	0.09
	7	0.19	0.21	0.19	0.09	0.09
	8	0.19	0.21	0.19	0.09	0.09
	9	0.19	0.21	0.19	0.09	0.09
	10	0.19	0.21	0.19	0.09	0.09
	11	0.19	0.21	0.19	0.09	0.09
	12	0.19	0.21	0.19	0.09	0.09

According to[START_REF] Claus | A plumbline constraint for the rational function lens distortion model[END_REF], the rational function model can also be solved linearly by minimizing some algebraic error. But this error is not directly related to the geometric error we want to minimize and sometimes leads to undesirable result.


Mediamobile Guillaume Allain

Fabrice Gamboa

Mediamobile Philippe Goudal

Jean-Noël Kien
email: jean-noel.kien@v-trafic.com

Modeling Weather Conditions Consequences on Road Trafficking Behaviors

Keywords: Forecasting method, Linear thresholded model, Road traffic, Spatial extrapolation, Weather

 to road trafficking. This model adapts itself locally to the whole road network and provides accurate unbiased forecasted speed using live or short term forecasted weather data information.

I. INTRODUCTION

I T is commonly accepted that adverse weather conditions modify significantly traffic flow dynamics in a complex way. Actually, it is well known that bad weather conditions such as, heavy rain, fog, snow, induce a significant decrease on traffic flow speeds. Note that it can be partially explained by the legal speed regulations. However, if several studies conclude that road traffic speed decreases during adverse weather, this trend is only confirmed. Furthermore, up to our knowledge, no quantitative analysis has been conducted to forecast the evolution of the observed speed of vehicles. Even deterministic models for road traffic fail in this study since they usually involve a large set of parameters which are all affected by the change of weather conditions. This prevents the use of such equations. In this paper, we tackle this issue and provide a general model to estimate the change in traffic flow speed for different weather scenarios.

Actually, trying to understand the impact of weather conditions requires a direct comparison between observed traffic speeds whose variability is only due to such changes. Hence, quantifying the impact of adverse weather conditions on traffic speed can only be conducted through the analysis of a two paired data. Each pair corresponds to similar traffic conditions but with different weather conditions. This is a difficult issue since road trafficking is a non stationary phenomenon and much of the variability is due to the road condition changes (with the occurrence of traffic jams). So, using empirical studies is not obvious because of the heterogeneity of the data. Moreover, some conditions are scarcely observed. This increases the difficulty of the estimation and requires to choose a set of weather conditions that makes sense for the road manager but also with a large enough number of observations. Some work has already been conducted in this direction. We refer to [START_REF] Faouzi | Toward Weather-Responsive Traffic Management on French Motorways[END_REF], [START_REF] El Faouzi | Motorway travel time prediction based on toll data and weather effect integration[END_REF], [START_REF] Faouzi | Effects of adverse weather on traffic and safety: State-of-the-art and a European initiative[END_REF], [START_REF] Billot | Multilevel Assessment of the Impact of Rain on Drivers' Behavior[END_REF], [START_REF] Billot | Integrating the impact of rain into traffic management : online traffic state estimation using sequential Monte Carlo techniques[END_REF] and [START_REF] Kilpeläinen | Effects of weather and weather forecasts on driver behaviour[END_REF]. Nevertheless, the drawback of previous methods is that the speed modifications are global. This means that the speeds are affected without enabling these changes to depend on the values of the initial velocity. In this work, we use a more flexible model using Multivariate Adaptive Regression Splines (MARS) model to consider thid initial condition. Furthermore, using a threshold enables to consider different impacts according to different levels of weather changes. Such procedure has been described in [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] and its implementation is detailed in [START_REF] Friedman | Fast MARS[END_REF]. Some work was conducted in the same direction, see for instance [START_REF] Garthwaite | Predicting multivariate responses in multiple linear regression -Discussion[END_REF], [START_REF] Mammen | Locally Adaptive Regression Splines[END_REF] and references therein. The calibration of the model is achieved by estimating the parameters on a learning set. This last set is built by selecting pairs of observations under the same traffic condition but with a different weather condition. We study the performance of this model and prove that it enables to forecast accurately the possible speed evolutions.

The paper falls into the following parts. Section II is devoted to the description of the data and their particularity. The following section, Section III describes the construction of the model while its performances are analyzed in Section IV. Finally, we discuss some of our results and draw some guidelines for speed forecast under adverse weather conditions in Section V.

II. DATA AND ISSUES

A. Description

A road network can be represented as a set of links connected together in a form depending on the underlying road network. Usually, links are classified by a well-known attribute: the Functional Road Class (FRC) [START_REF] Poovadan | MultiNet EUR 2011.12 Product Release Notes[END_REF]. FRC is a classification based on the importance of the road in the connectivity of the total network. Table I recaps the relation between the FRC and the network. In our study, industrial constraints make us work only on roads from FRC 0 to 2 because the most part of the information provided by Mediamobile that matches the customer demand concerns this FRC range. Vehicles equipped with some GPS device can return at each time their positions to a server. So that, they may be considered as floating sensors on the network. Such sensors form the Floating Car Data (FCD) source of traffic information. A map-matching algorithm (out of the scope of this article) establishes speed V(x, t) on a link x at a time t from a couple of successive positions and times by matching them on a digitized network. Then, when a significant number of speeds on a link x are at hand, we may produce a microscopic processed FCD speed V(x, t). By using this gathering technology, we are not geographically dependent of any counting station but we are limited by the GPS users feedback. Nevertheless, this source of traffic data is relevant since we have a huge amount of data (515 798 606 positions in March 2011 for instance) and we potentially cover all the network. The traffic database used in the paper was provided by Mediamobile and is composed of the vehicle speeds over T = 107 days.

The weather database used in this paper was also provided by Mediamobile through his partnership with Météo-France. Since 2009, Météo-France provides to Mediamobile a new high quality flow of data. It incorporates geo-tagged point every unit of 120 000 road kilometers of the France network. This flow of meteorological data is an aggregation of real time and forecasted observations. Then for our study, we have at hand a regular flow of weather data M(x, t) on a link x at a time t updated every 15 minutes. In this paper, we will focus on the following bad weather conditions: soft rain, medium and strong rain, rain and snow mixed, drizzle.

B. Data quality

Individual microscopic traffic data usually exhibit high noise and outliers due to several causes

• GPS logger accuracy, • incorrect vehicle path estimation: a wrong projection of the vehicle path can lead to incorrect speeds and further, speeds on incorrect links, • inner variations of individual vehicles in the traffic flow. To decrease the noise and eliminate outliers, we use a two step filtering algorithm:

1) the first step filters out aberrant data. Mediamobile estimates a Free Flow Speed (FFS) defined as the most likely speed in free flow traffic conditions. Then, we filter out aberrant data where speed records are higher than 150% of the reference speed i.e. the FFS, 2) the second step put out links that do not have enough records. Here, we arbitrarily fix to 100 measures the minimum amount of data necessary to keep a link.

After performing this algorithm, we have confident traffic data. Weather data are already consistent because they have been preprocessed by Météo-France. So that, they do not need any more treatment.

C. Issues

The main issues are twofold:

• building a learning set for weather condition. Actually, the main goal is to build couples of speeds at a given location observed under the same traffic condition but with a different weather condition in order to understand the weather conditions consequences on road trafficking behaviors. First, we need to associate both traffic data and weather data. The frequency of weather data flow is 15 minutes. So even if a weather condition is observed at a time t 0 only, we will propagate it to the whole interval [t 0 , t 0 + 15[ such that all speeds V(x, t) observed in this interval get paired with M(x, t 0 ), • finding a predictive rule to forecast velocities. We used the heuristic idea that adverse weather conditions do not affect the velocity in the same way. Indeed, we build a model that includes a different treatment for different ranges of speeds. This rule must be stable to be extended to the whole road network, yet providing good enough estimations.

III. MODELING

A. Regression

Our aim is to rectify the forecasted speed according to road weather conditions. Many methods exist to forecast speeds on a road network but this feature is out of the scope of this article. Some examples of such methods can be find in [START_REF] Loubes | Road trafficking description and short term travel time forecasting, with a classification method[END_REF] and [START_REF] Allain | A Statistical Framework for Road Traffic Prediction[END_REF]. In this work, we already have forecasted speed at hand and we want to correct them according to road weather conditions. This will be done by applying a correction on the forecasted speed. We focus on a bias depending on the speed. Indeed, the expertise of road trafficking theory shows that under adverse road weather conditions, drivers reduce their velocities when they go fast whereas they do not change it otherwise. An usual way to rectify the speed is obtained by using a polynomial function of degree 1

∀t ∈ V(t 0 ), V (x, t) = p k=1 P θ k V (x, t 0 ) × 1 [α k ,α k+1 ] V (x, t 0 ) (M 1 ) with θ k = [a k ; b k ], P θ k a polynomial of degree 1 such that P θ k V (x, t 0 ) = a k .V (x, t 0 ) + b k and level of speeds ∀k ∈ [1..p],
[α k ; α k+1 ] are such that α 1 = 0 and α p+1 = +∞. V (x, t 0 ) is the forecasted speed at time t 0 and V (x, t) its corrected speed. V(t 0 ) is a short term neighborhood around t 0 . This means that we will adjust V (x, t 0 ) in a time neighborhood while no newer speed data is available. This model is no more than the so-called Multivariate Adaptive Regression Splines (MARS) introduced in [START_REF] Friedman | Fast MARS[END_REF] where non linearities match driver's behavior changes. We refer to this model M 1 as the MARS model. Indeed, we can rewrite the model M 1 in a classical MARS form as

∀t ∈ V(t 0 ), V (x, t) = p k=1 Q Γ k h 1k , h 2k , h 3k with Γ k = [γ 1k ; γ 2k ; γ 3k ], Q Γ k a polynomial of degree 1 such that Q Γ k h 1k , h 2k , h 3k = γ 0 + γ 1k .h 1k + γ 1k .h 2k + γ 1k .h 3k , a set of hinge functions    h 1k = max(0, V (x, t 0 ) -α k ) h 2k = max(0, V (x, t 0 ) -α k+1 ) h 3k = max(0, α k+1 -V (x, t 0 )) and level of speeds ∀k ∈ [1..p], [α k ; α k+1 ] such that α 1 = 0 and α p+1 = +∞.
We point out in Section IV-B that the model M 1 suffers from a lack of stability on the network. Indeed, it has not an homogeneous structure among the network. This means that we are not able to extrapolate the model to all links. Moreover this model may apply a correction on low speeds under adverse weather condition. Hence, it does not fit the important feature: driving at low speed under bad weather conditions has no impact on drivers' behavior. That is why, we rather use the following variant of model M 1 that fits better the traffic flow theory and that is easier to extrapolate

∀t ∈ V(t 0 ), V (x, t) = V (x, t 0 ) × 1 [0, θ 0 1-θ 1 [ V (x, t 0 ) + θ 1 .V (x, t 0 ) + θ 0 × 1 [ θ 0 1-θ 1 ,+∞[ V (x, t 0 ) (M 2 )
We refer to model M 2 as the linear thresholded model since bad weather conditions have no consequence on road traffic flow under the critical speed θ 0 /(1 -θ 1 ). Beyond this value, vehicles are decreasing their speeds linearly.

B. An association between speeds

We wish to build a matching between speeds of vehicles at a given location observed under the same traffic condition but with a different weather condition. So that, we use the following scheme (illustrated in Fig. 1):

1) extract occurrence times of climate change t 0 ,

2) build a time neighborhood η(t 0 ) around t 0 in such a way that speeds in this neighborhood are stationary. In practice, we used to fix arbitrarily η(t 0 ) = [t 0 -h, t 0 ] with h = 5 minutes. This is generally narrow enough to warrant stationarity, 3) let t * be the time of the latest observed speed before the climate change. Finding at least one observation in η(t 0 ) for all t 0 is quite unlikely with our traffic data source. In fact, FCD are observed at random times so it is obvious that t * may not exist. The main consequence is a very sparse dataset. Thus, we decide to relax a little bit the assumption of traffic stationarity by allowing ourselves to pick similar velocities even if they belong to different observation days. This means that we assume the existence of a stationarity between days. The only criterion that matters to pair the data is their belonging to the same temporal neighborhood η(t 0 ). This ingenious consideration does not destroy more the stationarity than the climate change does itself, 4) associate V (x, t * ) to V (x, t 0 ) where M(x, t * ) = M(x, t 0 ). For instance, if we want to study the impact of rain on microscopic speeds, V (x, t * ) corresponds to speeds observed in no rain and no snow weather conditions (i.e. M(x, t * ) = NONE) and V (x, t 0 ) correspond to speeds observed in the rain (i.e. M(x, t 0 ) = RAIN).

Time Speed 0 20 60 100 t 0 η(t 0 )

Fig. 1. Matching between V (x, t * ) and V (x, t 0 ): the solid black line at t 0 corresponds to an occurrence of rain and V (x, t 0 ), the ringed star, is the corresponding speed; white points correspond to points included in η(t 0 ) = [t 0 -h, t 0 ] and V (x, t * ), the black point, corresponds to the latest speed in η(t 0 )

IV. RESULTS

We have introduced in Section III-A two models that are able to represent the impact of bad weather conditions on microscopic speeds. In this section, we first aim to select the best link by link model using a statistical approach. This gives a solution to a local modeling of our problem. Nevertheless, we also need to find a way to turn our link by link models into a network-wide model. This task is done right after the residual study for model selection.

The following results have been established with data presented in Section II matching the region of Toulouse in France during 107 days from November 1 st 2009 to February 15 th 2010. We also focus on the impact of the rain because it is the most common adverse weather condition in France. So we basically work with 256 832 observations on 2070 links.

A. Error

We highlighted two models: the MARS model M 1 and the linear thresholded model M 2 . Here, we are facing a classical problem in statistics. We have to choose only one of these two models. What is the best model? How can we compare models in order to select the best one? Such questions are classical in model selection area and a solution can be set up. We sample our data in two parts: a learning sample containing 90% of observations to estimate our models and a test sample containing the remaining 10% to validate and select the right model.

For each link of the network, we estimate both models by minimizing the Root Mean Squared Error (RMSE). The quality of a model is established by a classical goodness of fit value over all links calculated with the test sample

RM SE = x∈{links} 1 n x nx i=1 ( V (x, t) i -V (x, t) i ) 2
We obtain a RMSE of 10 359.8 for the MARS model and 11 795.27 for the linear thresholded model. So the MARS model fits the data better than the linear thresholded one but the difference of 1435.47 (13.8%) is not significant since it has been calculated on 552 links. Moreover, the MARS model has a more complex structure than the linear thresholded one because the second is nothing more than a particular constrained MARS model. So it fits the data better by construction. Nevertheless, RMSE's on each link are similar for both models as shown in Fig. 2. This means that although we conclude that MARS model is better, the linear thresholded model is not so far behind and have the undeniable advantage of being extrapolatable to all a network whereas MARS models cannot because they have not a homogeneous structure among the network. This will be detailed in the next section.

B. Stability and extrapolation

We aim to set up a simple method to apply a correction to microscopic speeds based on adverse weather conditions on all a complex network at a time rather than on each link of a network taken separately. In fact, the France road network from FRC 0 to 3 is composed of 1 740 462 links and obviously there is not the same number of different behaviors in response of adverse weather conditions. Thus, our will of finding a global method is justified. In this section, we discuss about how we can generalized link dependent models to a global network model. Interests of extrapolating our models match our industrial constraints:

• to sum up and simplify our link by link models, • to apply a correction on speeds under adverse weather conditions on uncovered links. First, we focus on the MARS model. For each link of the network, we have estimated a MARS model. We quickly face a problem to extrapolate this kind of model. Although this kind of model fits well the data, the structure is complex and different among links. For instance, Fig. 3 shows three theoretical MARS models on three links. We observe that the number of slopes and the number of points to model non-linearities may differ from one link to another, making a global structure for a network difficult to build. Moreover, this kind of model is not consistent with road trafficking theory. As a matter of fact, under adverse road weather conditions drivers may reduce their velocities for high level speeds while their behavior is unaffected for low ones which cannot be always the case in general. To sum up, link by link MARS models cannot be extrapolated to all a network mostly because they have not an homogeneous structure by link.

Second, we deal with linear thresholded models. Until now, our linear thresholded models were local because we built one model per link. With this kind of model, we always get a pair of parameters for each link which means that the structure is homogeneous among the network. So it is possible to extrapolate link by link models. Fig. 4 shows a 2-dimensional kernel density estimation of the distribution of the two parameters for all these models. Two modes appear: one is associated to low θ 0 's ( θ 0 ≤ 55 km/h) and another to high ones (θ 0 > 55 km/h). In fact, Fig. 5 shows exactly that marginal distribution for θ 0 is highly related to the FRC. We remind that we wish to build a global modeling for a network in other words our goal is to construct a model able to be applied on all the network. That is why we need to catch this dependence between θ 0 and the FRC. A way to do that is to find a normalization u such that the marginal distribution for θ0 = u(θ 0 , x) does not depend of the FRC. The use of u(•, x) = •/FFS(x) will normalize θ 0 correctly. Indeed, FFS is highly correlated with the FRC (see Fig. 6) and this normalization appears to be the most relevant and significant (5% significant F-tests has been done). Parameters of the global model remain as the empirical mean of the parameters of the link by link models i.e. ( θ 0 , θ 1 ) = ( θ0 , θ 1 ). One could naturally expect that (θ 0 , θ 1 ) is more than FRC dependent but also depend on the climate zone. Comparing the RMSE of the global model on climate zones with the RMSE of models aggregated by climate zone, we concluded that the global model was better in each case which means that the discrimination by climate zone is not significant.

Finally, our two link by link models M 1 and M 2 have similar RMSE on our data. Nevertheless, the MARS model can not be extrapolated by construction and it does not respect road trafficking theory. Thu,s the correct model for our problem is the linear thresholded model M 2 because it fits the data as the MARS model but it also respects road trafficking theory and it can be easily extrapolated to the whole network. Constraining the MARS model makes it homogeneous among all the links and consequently extrapolable. The form of the global linear thresholded model for all the network remains

∀t ∈ V(t 0 ), V (x, t) = V (x, t 0 ) × 1 [0, θ0 .F F S(x) 1-θ 1 [ V (x, t 0 ) + θ 1 .V (x, t 0 ) + θ0 .F F S(x) × 1 [ θ0 .F F S(x) 1-θ 1 ,+∞[ V (x, t 0 ) (M 3 ) It is interesting to rewrite the model M 3 in a interpretable form for road trafficking experts. Let α = θ0 1-θ1 and β = 1 -θ 1 , we obtain If V (x, t 0 ) ≥ α.F F S(x), V (x, t) = V (x, t 0 ) -β. V (x, t 0 ) -α.F F S(x) (M 4 ) Else V (x, t) = V (x, t 0 )
Thus the interpretation is really natural and illustrated in Fig. 8. If there is a vehicle at a speed above a proportion of the Free Flow Speed α.F F S(x) and it starts raining, we decrease its speed by β. V (x, t 0 ) -α.F F S(x) i.e. a proportion of the difference between its initial speed and the proportion of the Free Flow Speed . α.F F S(x) represents the speed at which adverse weather conditions start impacting speeds. We have built a global model which adapts itself locally since we use this ingenious normalization θ0 = θ 0 /F F S(x). With our data, estimates of M 3 are ( θ0 , θ 1 ) = (0.66, 0.16).

Let us consider a basic example to practice with the model: a vehicle is recorded at 130 km/h on a freeway where the Free Flow Speed is 130 km/h and it starts raining. So since 130 ≥ 0.66 × 130/(1 -0.16), we apply a correction and the speed is reduced to 130 -(016 × 130 + 0.66 × 130) = 106.6 km/h. Remark that it respects the French road speed limit on a freeway (130 km/h in general and decreasing to 110 km/h when raining).

Until now, we have selected the linear thresholded model and extrapolated it on a network. We know that our extrapolation by normalizing θ 0 by the Free Flow Speed was the best relevant according to our result but we also need to measure the actual loss of quality in extrapolating the model. So we calculate and compare these two following quantities still calculated on the test sample

RM SE |M2 = x∈{links} 1 n x nx i=1 ( V (x, t) i |M2 -V (x, t) i ) 2 RM SE |M3 = x∈{links} 1 n x nx i=1 ( V (x, t) i |M3 -V (x, t) i ) 2
with V (x, t) i |M2 the forecasted speed with a link by link linear thresholded model and V (x, t) i |M3 the forecasted speed with a linear thresholded model for the network i.e. link by link models extrapolated to all the network. We obtain RM SE |M2 = 11 795.27 and RM SE |M3 = 12 511.37. The loss of quality associated to the globalization ∆RM SE is equal to 6.07%. Thus we can warrant that the extrapolation of our link by link models is relevant and do not destroy the quality of fit compared to our local models.

V. CONCLUSION

We have tackled the issue of building a generic rule able to predict the evolution of vehicle velocities when the weather condition changes. For this, we have considered a version of Multivariate Adaptive Regression Spline model, well calibrated to get stable and accurate predictions. So we get model M 4 that basically does not modify speeds under a proportion of the free flow speed. Above this threshold value, the model decreases the speed by a proportion of the difference between this speed and the threshold value. Thus, one can use the model M 4 to correct forecasted speeds with the information of weather conditions. To learn the model over a data set, we had to construct a well adapted learning set. One of the difficulties of this task was to overcome the non stationarity of the observed data which mix both the variability due to the changes of the weather conditions and the one due to the changes in the traffic conditions. This was achieved by considering time neighborhood around similar velocities. Moreover the desired stability of the decision rule enables us to extend this method over a whole road network. This is, to our knowledge, the first global quantitative analysis of the impact of adverse weather on the observed vehicle velocities. This is a major improvement to forecast travel time with some knowledge on the weather conditions. This contributes to a better quality of the forecasts done by Mediamobile company.

Moreover, this study is also a key to better understand the macroscopic impact of adverse weather conditions on the Free Flow Speed. Actually, some fancy results emerge with the extrapolation of the linear thresholded model. We could think that adverse weather conditions not only impact microscopic speeds but also impact the well known Free Flow Speed. The global model M 3 includes such a result. When building this model, we build a speed at which adverse weather conditions start impacting speeds and this speed is nothing more than a proportion of the Free Flow Speed α.F F S(x). In this work, we provide another reference speed which can be considered as the Free Flow Speed under a certain adverse weather condition.

Fig. 2 .

 2 Fig. 2. Biplot of RMSE by links for MARS model and linear thresholded model

Fig. 3 .

 3 Fig. 3. Example of three MARS models on three links

Fig. 4 .Fig. 5 .

 45 Fig. 4. Global distribution of (θ 0 , θ 1 )

Fig. 6 .Fig. 7 .

 67 Fig. 6. Boxplots of Free Flow Speeds by FRC

Fig. 8 .

 8 Fig. 8. Percentage of speed reduction based on M 4

ACKNOWLEDGMENT

The authors would like to thank Météo-France and Mediamobile for providing respectively weather and road traffic data through their Road&Weather partnership.

This work was supported by Mediamobile with a partnership with Météo-France.


Julien Roth
email: julien.roth@univ-mlv.fr

A Remark on Almost Umbilical Hypersurfaces

Keywords: Hypersurfaces, Rigidity, Pinching, Ricci Curvature, Umbilicity Tensor, Higher Order Mean Curvatures. Mathematical Subject Clasification: 53A07, 53A10, 53C20, 53C24

In this article, we prove new stability results for almost-Einstein hypersurfaces of the Euclidean space, based on previous eigenvalue pinching results. Then, we deduce some comparable results for almost umbilical hypersurfaces.

Introduction

It is a well-known fact that a totally umbilical hypersurface of the Euclidean space R n+1 which is not totally geodesic is a round sphere. An other classical rigidty result states that an Einstein (with positive scalar curvature) hypersurface of the Euclidean space R n+1 is a round sphere. This was proved by Thomas [START_REF] Thomas | On closed space of constant mean curvature[END_REF] and independently by Fialkow [START_REF] Fialkow | Hypersurfaces of a space of constant curvature[END_REF] in the 30's. Recently, Grosjean [START_REF] Grosjean | Upper bounds for the first eigenvalue of the Laplacian on compact manifolds[END_REF] gave a new proof based on the equality case of an estimate of the first eigenvalue of the Laplacian involving the scalar curvature. It is then natural to consider the case of almost umbilical and almost Einstein hypersurfaces and ask if they are close to round spheres (in a sense to be precised). Shiohama and Xu [START_REF] Shiohama | Rigidity and sphere theorems for submanifolds[END_REF][START_REF] Shiohama | Rigidity and sphere theorems for submanifolds II[END_REF] proved that under a condition on Betti numbers, almost umbilical hypersurfaces of Euclidean space are homeomorphic to the sphere. For almost Einstein hypersurfaces, Vlachos obtained a comparable result in [START_REF] Vlachos | Almost-einstein hypersurfaces in the euclidean space[END_REF]. By studying the stablity problem of Grosjean's proof, we showed in [START_REF] Roth | Pinching of the first eigenvalue of the Laplacian and almost-Einstein hypersurfaces of Euclidean space[END_REF] another promity result for almost Einstein, namely these hypersurfaces are diffeomorphic and quasi-isometric to round sphere. We got stronger proximity but in counterpart, with stronger assumption that in the result of Vlachos. The aim of the present paper is to relax the assumptions of this result in one hand (Theorem 1) and on the other hand to get stability results for almost umbilical hypersurfaces to be compared with the results of Shiohama and Xu (see Theorems 2 and 3).

Preliminaries

Let (M n , g) be a n-dimensional compact, connected, oriented Riemannian manifold without boundary, isometrically immersed into the (n+1)-dimensional Euclidean space (R n+1 , can). The second fundamental form B of the immersion is the bilinear symmetric form defined by

B(Y, Z) = -g ∇ Y ν, Z ,
where ∇ is the Riemannian connection on R n+1 and ν the outward normal unit vector field on M .

From B, we can define the mean curvature,

H = 1 n tr (B),
and, more generally, the higher order mean curvatures,

H r = 1 n r σ r (κ 1 , • • • , κ n ),
where σ r is the r-th symmetric polynomial and κ 1 , • • • , κ n are the principal curvatures of the immersion. By convention, we set H 0 = 1. Note that H 1 = H and from the Gauss equation, H 2 is, up to a mutliplicative constant, the scalar curvature. Namely, we have H 2 = 1 n(n-1) Scal . These extrinsic curvatures satisfy the well-known Hsiung-Minkowski formula, for 1 r n,

M H r-1 + H r X, ν dv g = 0, (1)
where X is the position vector and ν the normal vector of the immersion. They also satisfy the following inequalities if H r is a positive function:

(2)

H 1 r r ≤ H 1 r-1 r-1 ≤ • • • ≤ H 1 2 2 ≤ H.
Moreover, Reilly [START_REF] Reilly | On the first eigenvalue of the Laplacian for compact submanifolds of Euclidean space[END_REF] proved some upper bounds for the first eigenvalue of the Laplacian for hypersurfaces of R n+1 in terms of higher order mean curvatures. Precisely, he showed

(3) λ 1 (M ) M H r-1 dv g 2 n Vol(M ) M H 2 r dv g ,
with equality if and only if m is a geodesic hyperspheres of R n+1 . By the Hölder inequality, we obtain for p 2,

λ 1 (M ) n ||H r || 2 2p Vol(M ) 2 M H r-1 dv g 2 .
Now, for p 2 and 1 r n, we define

k p,r = ||H r || 2 2p Vol(M ) 2 M H r-1 dv g
2 , which are the constants involved in Theorem 1.

Note that we use the following convention, for any smooth fonction f and any p 1,

||f || p = 1 Vol(M ) 1 p M |f | p dv g 1 p
.

The main tool that we will use in the present paper is the following pinching result, associated with these inequalities, that we proved in [START_REF] Roth | Pinching of the first eigenvalue of the Laplacian and almost-Einstein hypersurfaces of Euclidean space[END_REF].

Theorem A (Roth [START_REF] Roth | Pinching of the first eigenvalue of the Laplacian and almost-Einstein hypersurfaces of Euclidean space[END_REF]). Let (M n , g) be a compact, connected, oriented Riemannian manifold without boundary isometrically immersed in R n+1 , n 2. Let r ∈ {1, • • • , n} and assume that H r > 0 if r > 1. Then for any q 2 and any θ ∈]0, 1[, there exists a constant K θ depending only on n, ||H|| ∞ , Vol(M ), ||H r || 2q and θ such that if the pinching condition

(P K θ ) 0 λ 1 (M ) M H r-1 dv g 2 -nVol(M ) 2 ||H r || 2 2q > -K θ is satisfied, then M is diffeomorphic and θ-quasi-isometric to S n n λ 1 .
Remark 1. Note that if q n 2r ; then K θ does not depend on ||H r || 2q . Remark 2. The case r = 1 was proved by Colbois and Grosjean [START_REF] Colbois | A pinching theorem for the first eigenvalue of the Laplacian on hypersurfaces of the Euclidean space[END_REF]. In that case we do not need to assume that H > 0.

Almost Einstein hypersurfaces

We showed in [START_REF] Roth | Pinching of the first eigenvalue of the Laplacian and almost-Einstein hypersurfaces of Euclidean space[END_REF] that almost-Einstein hypersurfaces of R n+1 are close to round spheres. Namely, Theorem B (Roth [START_REF] Roth | Pinching of the first eigenvalue of the Laplacian and almost-Einstein hypersurfaces of Euclidean space[END_REF]). Let (M n , g) be a compact, connected, oriented Riemannian manifold without boundary isometrically immersed in R n+1 , n 2.

Let θ ∈]0, 1[. If (M n , g) is almost-Einstein, that is, ||Ric -(n -1)kg|| ∞ ε for a positive constant k, with ε small enough depending on n, k, ||H|| ∞ and θ, then M is diffeomorphic and θ-quasi-isometric to S n 1 k
By θ-quasi-isometric, we understand that there exists a diffeomorphism

F from M into S n 1 k
such that, for any x ∈ M and for any unitary vector u ∈ T x M , we have

|d x F (u)| 2 -1 θ.
This theorem is a corollary of our pinching result for the first eigenvalue of the Laplacian (Theorem A).

In this article, we consider almost-Einstein hypersurfaces of R n+1 in a weaker sense, namely for the L q -norm, that is, ||Ric -(n -1)kg|| q ε for some positive constant k and a sufficiently small ε. We prove that for some suitable constants k, such manifolds are close to round spheres. Precisely, we prove the following Theorem 1. Let (M n , g) be a compact, connected, oriented Riemannian manifold without boundary isometrically immersed in

R n+1 . Let q > n 2 , r ∈ {1, • • • , n} and if r > 1, assume that H r > 0. Let θ ∈]0, 1[, if (M n , g) satisfies ||Ric -(n -1)k p,
r g|| q ε for some sufficiently small ε depending on n, q, ||H|| ∞ , Vol(M ) and θ, then M is diffeomorphic and θ-quasi-isometric to S n 1 kp,r .

The constants k p,r in the theorem are defined from the higher order mean curvature (see Sect. 2).

After giving the proof of this theorem, we will give some applications to almost-umbilical hypersurfaces.

Proof: The proof is based on the above pinching result combined with a lower bound for the first eigenvalue of the Laplacian due to Aubry [START_REF] Aubry | Finiteness of π 1 and geometric inequalities in almost positive Ricci curvature[END_REF]. Assume that ||Ric -(n -1)kg|| q ε(n, q, k) for a positive constant k, q > n 2 and ε small enough, then from Theorem 1.1 of [START_REF] Aubry | Finiteness of π 1 and geometric inequalities in almost positive Ricci curvature[END_REF], we deduce that λ 1 (∆) satisfies

(4) λ 1 (∆) nk(1 -C ε ),
where C ε is an explicit constant such that C ε -→ 0 when ε -→ 0. Now, with the particular choice of k = k p,r , we get:

λ 1 (M ) M H r-1 dv g 2 -nVol(M ) 2 ||H r || 2 2p > -K ε for some constant K ε such that K ε -→ 0 when ε -→ 0.
Let θ ∈]0, 1[, we choose ε(n, q, k, θ) small enough such that K ε is small enough in Theorem A to obtain a diffeomorphism and θ-quasi-isometry between M and S n 1 kp,r .

Remark 3. Note that in this Theorem, ε depends on Vol(M ) contrary to Theorem B. We can remove this dependence on the volume and replace it by a dependence on ||H r-1 || 1 . Indeed, using an upper bound on the volume under the L p condition on Ricci (result by Aubry [START_REF] Aubry | Finiteness of π 1 and geometric inequalities in almost positive Ricci curvature[END_REF]) we have that the volume is bounded from above by a constant depending on n, q, ||H|| ∞ and k r,p that is in fact on n, q, ||H|| ∞ and ||H r-1 || 1 because of the definition of k r,p and (2). On the other hand, by the classical extrisnic Sobolev inequality of Hoffman and Spruck [START_REF] Hoffman | Sobolev and isoperimetric inequalities for Riemannian submanifolds[END_REF], we can bound the volume from below by a constant depending only on n and ||H|| ∞ . Now, we will deduce from Theorem 1 some Corollaries for almost-umbilical hypersurfaces of R n+1 .

Almost umbilical hypersurfaces

First, we give the following theorem, which a direct application of Theorem B.

Theorem 2. Let (M n , g) be a compact, connected, oriented Riemannian manifold without boundary isometrically immersed in

R n+1 . Let θ ∈]0, 1[. If (M n , g) is almost-umbilical, that is, ||B -kg|| ∞ ε
for a positive constant k, with ε small enough depending on n, k and θ then M is diffeomorphic and θ-quasi-isometric to S n 1 k . Proof : From the Gauss formula, we have

(5) Ric(Y, Y ) = nH B(Y ), Y -B(Y ), B(Y ) ,
for a tangent vector field Y . From (5) and ||B -kg|| ∞ ε, we deduce

Ric(Y, Y ) nk 2 ||Y || 2 (1 -ε) 2 -k 2 ||Y || 2 (1 + ε) 2 (n -1)k 2 ||Y || 2 -α n (ε)||Y || 2 ,
where α n is a positive function such that α n (ε) -→ 0 when ε -→ 0. Similarly, we get

Ric(Y, Y ) (n -1)k 2 ||Y || 2 + α n (ε)||Y || 2 .
Finally, we have ||Ric -

(n -1)k 2 g|| ∞ α n (ε),
which implies, by Theorem B, that for ε small enough, M is diffeomorphic and θ-quasi-isometric to S n 1 k . Now, from Theorem 1, it is possible to obtain, in some particular cases, comparable results for almost umbilical hypersurfaces in an L q -sense. We recall that the umbilicity tensor is defined by τ = B -HId .

As we mentionned above, if M is umbilical, i.e.,τ = 0, and if M is compact, then M is a geodesic sphere. Here, we prove the following stability result for almost umbilical hypersurfaces.

Theorem 3. Let (M n , g) be a compact, connected, oriented Riemannian manifold without boundary isometrically immersed in R n+1 . Let q > n 2 , r ∈ {1, • • • , n} and if r > 1, assume that H r > 0. For any θ ∈]0, 1[, there exists two constants

ε i (θ, n, ||H|| ∞ , Vol(M )), i = 1, 2, such that if 1. ||τ || 2q ε 1 , 2. ||H 2 -k p,r || q ε 2 , for p 4 and 1 r n, then M is diffeomorphic and θ-quasi-isometric to S n 1 √ kp,r
Remark 4. Note that for r = 1, the result is due to Grosjean-Roth (see [START_REF] Grosjean | Eigenvalue pinching and application to the stability and the almost umbilicity of hypersurfaces[END_REF]), using a pinching result which involves only H 1 (see [START_REF] Colbois | A pinching theorem for the first eigenvalue of the Laplacian on hypersurfaces of the Euclidean space[END_REF]).

Proof. Here again, from the Gauss formula, we have

Ric = nHB -B 2 .
From this, we deduce that

Ric -(n -1)H 2 g = nHB -B 2 -(n -1)H 2 g = (n -2)Hτ -τ 2 , which implies ||Ric -(n -1)kg|| q ||Ric -(n -1)H 2 g|| q + (n -1) √ n||(H 2 -k)|| q (n -2)||H|| ∞ ||τ || 2q + ||τ || 2 2q + (n -1) √ n||(H 2 -k)|| q (n -2)||H|| ∞ ε 1 + ε 2 1 + (n -1) √ nε 2
Now, we conclude by taking ε 1 and ε 2 small enough depending on n, ||H|| ∞ and θ in order to apply Theorem 1 and obtain the θ-quasi-isometry.

Then, we deduce the following corollary which is to compare with Theorem 2 . For ε small enough, the assumptions of Theorem 3 are satisfied and we can conclude that M is diffeomorphic and quasi-isometric to S n 1 √ kp,r . Remark 5. We want to point out that this corollary is an improvement of Theorem 2 only in some sense. Indeed, we improve the L ∞ -proximity to an L 2q -proximity, but this corollary is valid only for some special constants k p,r and not for any positive constant as in Theorem 2.

Corollary 1 . 1 √

 11 Let (M n , g) be a compact, connected, oriented Riemannian manifold without boundary isometrically immersed in R n+1 . Let θ ∈]0, 1[.If (M n , g) is almost-umbilic, that is, ||Bk p,r g|| 2qε, for q > n 2 , with ε small enough depending on n, ||H|| ∞ and θ then M is diffeomorphic and θ-quasi-isometric to S n kp,r . Proof : A simple computation shows that||H 2k p,r || 2q α 1 ||Bk p,r g|| 2q , and ||τ || 2q α 2 ||Bk p,r g|| 2q , for two constants α 1 and α 2 depending on n and ||H|| ∞ . Since we assume that ||Bk p,r g|| 2q ε, we get 1. ||H 2k p,r || 2q α 1 ε, 2. ||τ || 2q α 2 ε.