Maxime Lhuillier
email: maxime.lhuillier@univ-bpclermont.fr

Automatic Scene Structure and Camera Motion using a Catadioptric System

Keywords: Central and Non-Central Catadioptric Cameras, Structure from Motion, Bundle Adjustment, Vision System

Fully automatic methods are presented for the estimation of scene structure and camera motion from an image sequence acquired by a catadioptric system. The first contribution is the design of bundle adjustments for both central and non-central models, by taking care of the smoothness of the minimized error functions. The second contribution is an extensive experimental study for long sequences of catadioptric images in a context useful for applications: a hand-held and equiangular camera moving on the ground. An equiangular camera is non-central and provides uniform resolution in the image radial direction. Many experiments dealing with robustness, accuracy, uncertainty, comparisons between both central and non-central models, and piecewise planar 3D modeling are provided.

Introduction

The automatic estimation of scene structure and camera motion from image sequence acquired by a perspective camera (whether calibrated or not) taken by hand has been a very active topic [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Faugeras | The Geometry of Multiple Images[END_REF] during recent decades, and many successful systems now exist [START_REF]d3 Ltd[END_REF][START_REF] Pollefeys | Self-Calibration and Metric Reconstruction in spite of Varying and Unknown Internal Camera Parameters[END_REF][START_REF] Nister | Visual Odometry[END_REF][START_REF] Lhuillier | A Quasi-Dense Approach to Surface Reconstruction from Uncalibrated Images[END_REF]. This paper (an extended version of [START_REF] Lhuillier | Automatic Structure and Motion using a Catadioptric Camera[END_REF]) focuses on this estimation using a catadioptric system/camera. Such a system has a wide field of view thanks to a convex mirror mounted in front of a perspective camera as shown in Figure 1. A catadioptric system may be "central" if all back-projected rays go through a single point in 3D called the "center". An overview of the paper and comparisons with previous works are given in Sections 1.1 and 1.2, respectively. Section 1.3 provides a brief summary on the known desirable conditions for errors minimized by bundle adjustments and discusses errors introduced in previous works.

Overview

First, central and non-central catadioptric camera models are described in Section 2. The non-central model explicitly involves the ray reflection onto the mirror surface, while the central model directly defines a mapping from 3D to 2D. These models require neither a conic as mirror profile nor the perspective camera center at one of the conic focal points.

Second, Section 3 presents bundle adjustments for a catadioptric camera in a general context: field of view greater than 180 • , scene with distant 3D points, approximate calibration knowledge and image noise. Although bundle adjustment is a well known iterative method, we spend some time in Sections 3.1 and 3.2 going over desirable smoothness conditions on the minimized errors to facilitate the convergence. These conditions are not fully satisfied by the standard 2D reprojection error of a catadioptric camera. With this in mind, we propose bundle adjustment methods minimizing 2D image errors (measured in pixels) in Section 3.3 and angular errors (measured in radians) in Section 3.4.

Third, Section 4 describes the automatic estimation method. The geometry of the image sequence is estimated for the central model. Then, a second estimation using the non-central model is obtained. In both cases, an approximate and given calibration is refined with the 3D. Fourth, Section 5 presents an extensive experimental study in a context that we expect to be representative (although not trivial) and useful for applications: a hand-held and equiangular catadioptric camera moving on the ground. An equiangular camera is non-central and provides uniform resolution in the image radial direction. We propose experiments about robustness, accuracy and uncertainty estimation, performance comparisons between central and non-central models, and piecewise planar 3D modeling from the reconstructed points.

The two main contributions of this paper are the bundle adjustments with image and angle errors in Section 3, and the experimental study in Section 5.

Comparisons with Previous Works

The most related work is that of Micusik and Pajdla [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF]. It is a two step method like ours: (1) estimate the sequence geometry using a central model for the camera and (2) upgrade the sequence geometry using a non-central model enforcing the mirror knowledge. The central model is similar to ours and is more general than the Geyer-Daniilidis model [START_REF] Geyer | A unifying Theory for Central Panoramic Systems and Practical Implications[END_REF] which requires a conic as mirror profile and the perspective camera center at one of the conic focal points. A general central model is more adequate to approximate a (general) non-central catadioptric camera. The differences are the following. The noncentral model proposed in [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF] is formalized from 2D to 3D since the authors minimize a 3D error for bundle adjustment. Our 3D to 2D formalization is more adequate for bundle adjustment minimizing a 2D image error (measured in pixels). Comparisons between 3D and 2D errors are made in Section 1.3. Furthermore, the work [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF] emphasizes auto-calibration in a general context for two views only. This is not our focus: we assume that an approximate calibration is given and we emphasize [START_REF] Adiv | Inherent ambiguities in recovering 3-D motion and structure from a noisy flow field[END_REF] the definition of errors with high smoothness to be minimized by multi-view bundle adjustment and (2) extensive experiments on long sequences. From a technical viewpoint, we use the 7-point algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] to estimate fundamental matrix instead of the 9 and 15-point methods [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF] involving polynomial eigenvalue problems. Both papers provide experiments with a non-central catadioptric camera: a (roughly) orthographic camera with a parabolic mirror [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF], and a (roughly) equiangular catadioptric camera in our case. The mirror profile of our equiangular camera is not a conic.

A different problem [START_REF] Aliaga | Accurate Catadioptric Calibration for Real-time pose estimation in Room-size Environments[END_REF] is the real-time pose estimation of a parabolic catadioptric camera for planar motions in a room-size environment, using a few known 3D beacons. The ideal model is central: a parabolic mirror in front of an orthographic camera. The real model is non-central: a parabolic mirror in front of a perspective camera. In this context, the pose accuracy is improved by the non-central model. In our context, experiments show that the advantages of the non-central model are not so convincing.

The geometry of a catadioptric sequence without accurate calibration has been estimated in other works, but they involve two views and use central camera models. Calibration and successive essential matrices are estimated for a parabolic catadioptric camera from tracked points in the image sequence [START_REF] Kang | Catadioptric Self-Calibration[END_REF]. Calibration initialization is given by the approximate field of view angle. Other work on the same sensor [START_REF] Geyer | Structure and Motion from Uncalibrated Catadioptric Views[END_REF] introduces the parabolic fundamental matrix, which defines a bilinear constraint between two matched points represented in an adequate space, and shows that calibration and essential matrix recovery are possible from this matrix.

Which Errors for Bundle Adjustment ?

Bundle adjustment (BA) is the standard process to obtain an accurate estimate of the geometry (all cameras and 3D points) by minimizing a sum of squares of non-linear errors. A survey about bundle adjustment [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] (also refer to [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Madsen | Methods for Non-Linear Least Squares Problems Problems[END_REF]) describes the desiderata for these errors in detail.

First, errors should be smooth enough in order to facilitate the convergence of BA. C 2 continuity is recommended to reach quadratic final convergence using the Levenberg-Marquardt method involved in BA [START_REF] Madsen | Methods for Non-Linear Least Squares Problems Problems[END_REF][START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. This continuity includes the case of a 3D point at the infinity plane thanks to the homogeneous coordinates [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. It is easy to satisfy this condition for a conventional (perspective) camera, but it is not for a catadioptric camera. Indeed, the definition of errors with high smoothness is not straightforward in our context since the catadioptric projection function is always discontinuous at the infinite plane.

Second, errors should be well chosen in order to improve the geometry estimation in a statistical sense [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. Minimized errors are usually noisy measurements of physical quantities and are modeled by realizations of random vectors. We obtain a Maximum Likelihood Estimation (MLE) of the geometry with the common statistical model [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]: errors obey zero-mean isotropic Gaus-sian probability distributions which are independent and identical (outliers are removed before the BA step). The minimization of 2D reprojection errors has been used for many decades, whereas 3D error minimizations [START_REF] Kang | 3D Scene Data Recovery Using Omnidirectional Multibaseline Stereo[END_REF][START_REF] Lu | Fast and Globally Convergent Pose Estimation from Video Images[END_REF][START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF] provide biased results as mentioned in [START_REF] Kang | 3D Scene Data Recovery Using Omnidirectional Multibaseline Stereo[END_REF][START_REF] Lu | Fast and Globally Convergent Pose Estimation from Video Images[END_REF]. The reason is the following: the common Gaussian model is not tenable for 3D errors (distances between back-projected rays and 3D reconstructed points) because the true probability distribution of 3D error depends too much on the distance between cameras and the point. However, the bias is limited in experiments [START_REF] Kang | 3D Scene Data Recovery Using Omnidirectional Multibaseline Stereo[END_REF][START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF] since 3D points in indoor scenes have distances of the same magnitude order. More details on the 3D error in papers [START_REF] Kang | 3D Scene Data Recovery Using Omnidirectional Multibaseline Stereo[END_REF][START_REF] Lu | Fast and Globally Convergent Pose Estimation from Video Images[END_REF][START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF] are given below.

The 3D error is used for the point estimation by intersection of rays [START_REF] Kang | 3D Scene Data Recovery Using Omnidirectional Multibaseline Stereo[END_REF] (page 181). Let r(λ k ) = t k + λ k v k be the parametrization of the k th ray with t k and v k the starting point and normalized direction of the ray. The 3D point p closest to all of the rays minimizes k ||e k || 2 with the error

e k = p -(t k + λ k v k ).
No iterative method is needed since a closed form expression of p exists. The authors mention that they obtain a more "optimal" result by minimizing k

1 λ 2 k ||e k || 2 . They introduce the weights 1 λ 2 k
and argue that the uncertainty in point location grows linearly with λ k . In other words, the common Gaussian model above is more tenable for e k λ k than e k . Furthermore, we note that

|| e k λ k || = || p-t k λ k -v k ||
at the solution is similar to a distance between two close points onto the unit sphere: it is an angle between two directions. Section 3.4 describes an angle error which is C 2 continuous, even at the infinite plane. The errors of [START_REF] Kang | 3D Scene Data Recovery Using Omnidirectional Multibaseline Stereo[END_REF] are discontinuous at the infinite plane.

The 3D error is used for the estimation of the camera pose [START_REF] Lu | Fast and Globally Convergent Pose Estimation from Video Images[END_REF]: find the pose given many 3D points in the world coordinate system and the corresponding ray directions in the camera coordinate system. The method is globally convergent, but the authors mention and experiment that the pose solution more heavily weights the points that are farther away from the camera.

Last, the 3D error is used for bundle adjustment with a non-central catadioptric camera [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF]. In this context, the 3D error is preferred to the 2D reprojection error which is more time consuming. The reason is that the catadioptric projection function has no closed form for a general mirror. The bundle adjustment is designed for two views only. The cost function j e 2 j is minimized with error e j = (t 0 jt 1 j ).

v 0 j ×v 1 j ||v 0 j ×v 1 j || .
In this expression, the j th reconstructed point is in the center of the shortest transversal of two back-projected rays defined by (t 0 j , v 0 j ) and (t 1 j , v 1 j ) (t i j and v i j are the starting point and normalized direction of the back-projected ray by the i th camera). We see that the point parameters are eliminated. Drawbacks are the 2-view use, and the bias mentioned by previous authors for point [START_REF] Kang | 3D Scene Data Recovery Using Omnidirectional Multibaseline Stereo[END_REF] and camera pose [START_REF] Lu | Fast and Globally Convergent Pose Estimation from Video Images[END_REF] estimations.

Catadioptric Camera Models

Figure 5 shows a catadioptric camera we model: the image projection of the scene is between a large circle and a small circle. Sections 2.1 and 2.2 present the projection functions of the central and non-central camera models for a finite 3D point, respectively. Section 2.3 introduces the definition of antipodal projections. The case of a point in the infinite plane is detailed in Section 2.4.

Central Model

The central catadioptric model is defined by its orientation R (a rotation), the center t ∈ 3 ( is the real numbers), both expressed in the world coordinate system, and a projection function C from 3 to 2 defined below. If X is the homogeneous world coordinates of a finite 3D point such that the last coordinate is positive, the projection of X in the image is defined by

p c (X) = C(d), d = R [I 3 | -t]X. ( 1
)
Vector d is the direction of the ray from t to X in the camera coordinate system. The sign constraint for X is necessary to choose a ray (half line) among the opposite rays defined by directions d and -d. Indeed, these two rays are projected on two different points in the omnidirectional image if the field of view is greater than 180 • .

The central model also has a symmetry around the z-axis of the camera coordinate system: the omnidirectional image is between two concentric circles of radii r i up and r i down , and C is such that there is a positive and decreasing function r such that C(x, y, z) = r(α(x, y, z))

1 √ x 2 + y 2 x y . (2)
Let α = α(x, y, z) = α(d) be the angle between the z-axis and the ray direction d. The two angles α up and α down which define the field of view are such that α up ≤ α ≤ α down and r i down ≤ r(α) ≤ r i up (see Figure 2). We assume that r and p c are C 2 continuous functions.

Non-Central Model

The non-central catadioptric model is defined by the orientation R (a rotation) and the location t ∈ 3 of the mirror coordinate system, both expressed in the world coordinate system, the intrinsic parameters K p , orientation R p and location t p of the perspective camera expressed in the mirror coordinate system, and the (known) mirror surface.

The mirror surface has a symmetry around the z-axis of the mirror coordinate system. t is at the mirror apex. Let M + (A, B) be the function which gives the reflection point on the mirror for the ray which goes across points A and B (the mirror profile is such that there is at most one reflection point). Point A is in the incident part and B is in the reflected part of the ray, or vice versa. We have A, B, M + (A, B) ∈ 3 , all expressed in the mirror coordinate system. This function does not have a closed form, and its calculation from the mirror surface is presented in Appendix D.

Let X be the homogeneous world coordinates of a finite 3D point and

π 2 (( x y z ) ) = ( x z y z ) , π 3 (( x y z t ) ) = ( x t y t z t ) . (3) Since R (π 3 (X) -t) and K p R p [I 3 | -t p ]
are respectively the X coordinates and the perspective camera matrix in the mirror coordinate system, the projection of X in the image is defined by

p + nc (X) = π 2 (K p R p (M + (t p , R (π 3 (X) -t)) -t p )). (4)
Figure 3 shows a cross section of the mirror and pinhole camera, the points t p , X, p + nc (X), M + (t p , X) with R = I 3 and t = 0. Figure 3 also defines a kind of antipodal mirror point M -(t p , X) for M + (t p , X) and its projection p - nc (X) which will be useful later. Obviously, we have Fig. 3. The ray which goes across finite point X and pinhole camera center t p is reflected on the mirror at point M + (t p , X), and projected by the pinhole camera to p + nc (X). We also define the point M -(t p , X) on the mirror, projected by the pinhole camera to p - nc (X). The only difference between M -(t p , X) and M + (t p , X) is the following: the reflected ray from M + (t p , X) (respectively, M -(t p , X)) towards X is pointing outside (respectively, inside) the mirror. In general, X, M + (t p , X) and M -(t p , X) are not collinear points.

p - nc (X) = π 2 (K p R p (M -(t p , R (π 3 (X) -t)) -t p )). ( 5
)
We note that the parametrization (R p , t p , R) is not minimal because of the mirror symmetry around the z-axis: the expressions of p + nc (X) and p - nc (X) are unchanged by replacing (R p , t p , R) by (R z R p , R z t p , RR z ) and applying the mirror symmetry relations

R z M + (A, B) = M + (R z A, R z B), R z M -(A, B) = M -(R z A, R z B)
for any rotation R z around the z-axis. Last, we assume in this paper that the mirror profile is such that M + and M -are C 2 continuous functions. The resulting p + nc (X) and p - nc (X) are C 2 continuous.

Antipodal Projections

The antipodal projections of X are defined by {p c (X), p c (-X)} in the central case, and by {p + nc (X), p - nc (X)} in the non-central case.

Assume that both antipodal projections exist. In the central case, the antipodal projections are C(d) and C(-d) with d ∈ 3 . They are separated by the small circle. In the non-central case, Figure 3 shows that the small circle also separate the two antipodal projections. In practice, the distance between the two antipodal points is greater than the diameter of the small circle (hundreds of pixels in the experiments).

Catadioptric Cameras and the Infinite Plane

Let X ∞ = ( x y z 0 ) be a point in the infinite plane. We consider the finite point X(t) = ( x t y t z t 1 ) such that t ∈ converges to 0. The limit of X(t) is the infinite point X ∞ in the projective space, and two different limits of p c (X(t)) or p + nc (X(t)) are possible in the catadioptric image: one for each sign of t. These limits (if they exist) define the projection(s) of X ∞ by the catadioptric camera.

The main subject of this part is to show that the limit of each antipodal projection (defined in Section 2.3) of X(t) is a projection of X ∞ .

Obviously, the limits of the first antipodal projections

lim t → 0 t > 0 p c (X(t)), lim t → 0 t > 0 p + nc (X(t))
are projections of X ∞ for the central and non-central models, respectively.

In the central case, we introduce

d(t) = R [I 3 |-t]X(t) ||R [I 3 |-t]X(t)|| which verifies p c (X(t) = C(d(t)), p c (-X(t) = C(-d(t)), lim t → 0 t < 0 d(t) = -lim t → 0 t > 0 d(t).
Thanks to the continuity of C onto the unit sphere,

lim t → 0 t > 0 p c (-X(t)) = C(-lim t → 0 t > 0 d(t)) = C( lim t → 0 t < 0 d(t)) = lim t → 0 t < 0 p c (X(t)).
Thus, the limit of the second antipodal projection (on the left) is the second projection of X ∞ (on the right).

In the non-central case, a proof sketch is given. Figure 4 is obtained from Figure 3 for a finite point X = X(t) which converges to the right toward

X ∞ with t > 0. The points M + (t p , X), p + nc (X), M -(t p , X), p - nc (X) converge to M + (t p , X + ∞ ), p + nc (X + ∞ ), M -(t p , X + ∞ ), p - nc (X + ∞ ), respectively. We see in Fig- ure 4 that the back projected ray of pixel p - nc (X + ∞ ) is reflected by the mirror at M -(t p , X + ∞ )
with the direction given by X ∞ toward the left. Thus, the limit p - nc (X + ∞ ) of the second antipodal projection p - nc (X(t)) is one of the two projections of X ∞ . Furthermore, it is not p + nc (X + ∞ ), the first projection of X ∞ .

+ (t p , X + ∞ ) and M -(t p , X + ∞ ) have the same direction X + ∞ .
3 Bundle Adjustments for a Catadioptric Camera

Bundle adjustment [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] or BA is the standard process to obtain an accurate estimate of the sequence geometry. It iteratively adjusts the parameters of all cameras and all 3D points in order to minimize a cost function. This function is a sum of squares of non-linear errors

E = i,j ||e i j || 2 (6)
with e i j the error corresponding to the i th camera and the j th point. ||x|| is the Euclidean norm of a vector x. The standard 2D reprojection error is e i j = p i (X j )u i j (measured in pixels), where p i is the projection function by the i th camera, X j is the j th 3D point in the world coordinate system, and u i j the point detected in the i th image which corresponds to the j th point. Sections 3.1 and 3.2 describe smoothness conditions on e i j . Then, we introduce the image error in Section 3.3 and the angular error in Section 3.4.

Smoothness Conditions

The most commonly used minimization tool for geometry refinements is the Levenberg-Marquardt method (LM), which combines the Gauss-Newton algorithm and the descent of gradient [START_REF] Madsen | Methods for Non-Linear Least Squares Problems Problems[END_REF]. Thus, a summary of LM is useful to justify conditions on e i j . LM requires a starting point and produces a series of vectors which converges to a local minimizer of the cost function. This local minimizer is the right one if the starting point is close enough. Let p n be the current value of the geometry parameters. The quadratic Taylor expansion of the cost function defined by Equation 6is

E(p n + ∆) ≈ E(p n ) + ∆ E n + 1 2 ∆ E n ∆
with E n and E n the gradient and Hessian of E at p n . One LM iteration is

p n+1 = p n + ∆ n where ∆ n satisfies (E n + λ n I)∆ n = -E n
, and λ n is a well chosen parameter evolving during the convergence such that E(p n + ∆ n ) < E(p n ). I is the identity matrix (or a diagonal matrix obtained by the diagonal coefficients of E n [START_REF] Press | Numerical Recipes in C[END_REF]). If λ n is "small", the LM behavior is Newton like: we get quadratic convergence in the immediate neighborhood of the solution if E is C 2 continuous and E n is definite positive. Elsewhere, a gradient descent like behavior is preferred with a "large" λ n to guarantee a decreasing cost function.

In practice, the Gauss-Newton approximation is always used for the Hessian E n of least squares:

1 2 ∂ 2 (||e i j || 2 ) ∂p 0 ∂p 1 = e i j ∂ 2 e i j ∂p 0 ∂p 1 + ∂e i j ∂p 0 ∂e i j ∂p 1 ≈ ∂e i j ∂p 0 ∂e i j ∂p 1 (7)
where p 0 , p 1 are any parameters of the geometry. This approximation is acceptable if ||e i j || is small and ||

∂ 2 e i j ∂p 0 ∂p 1 || is bounded.
This summary about LM gives many conditions for an ideal error. First, the C 1 continuity of e i j is sufficient for the evaluation of e i j Jacobian (Equation 7). Second, the C 2 continuity of e i j is even better to reach the quadratic final convergence, with the additional condition: e i j = 0 if the i th camera and the j th point are in exact agreement. Third, the Jacobians of e i j should not all be 0 if e i j = 0. Otherwise, the E n approximation by Equation 7converges to zero at the final convergence of LM and the quadratic convergence is lost.

3D Point Parametrization and Crossing the Infinite Plane

In this part, we spend some time explaining why we should also pay attention to the parametrization of the 3D point X j involved in e i j in the neighborhood of the infinite plane. The case of e i j defined by the standard 2D reprojection error is studied (a similar discussion on this topic is also given in [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]).

Assume that the X j parametrization is (x, y, z). Points onto the infinite plane are ignored by this parametrization, although they may have a projection by the camera. Furthermore, major changes of x, y and z are needed to modify ||e i j || if the 3D point X j is distant from the i th camera. Now, the affine parametrization (x, y, z) is replaced by the homogeneous projective parametrization (x, y, z, t). This parametrization allows X j to go across (and onto) the infinite plane with small variations of t and x, y, z. After a crossing of the infinite plane (if the sign of t changes during an LM iteration), a point in front of the camera goes behind the camera and vice versa. A crossing of the infinite plane is just a kind of "faster way" to connect two distant 3D points of a camera. We note that e i j for a perspective camera is C 2 continuous with the 3D point parametrization (x, y, z, t), even if X j is in the infinite plane.

As mentioned in [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF], such a crossing of the infinite plane by a point is sometimes necessary for the bundle adjustment convergence due to noise or approximate calibration. An approximate calibration or noise may reconstruct many distant points (by ray intersection) behind the camera although they should be in front of it, and the only way to escape from this configuration using an LM iteration may be by crossing the infinite plane. Furthermore, such transitory configurations help with LM convergence. Now, we focus on the case of any catadioptric camera with e i j defined by the standard 2D reprojection error and the parametrization (x, y, z, t). This was not done in [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. A 3D point in the infinite plane is also visible, and it may have two projections by the camera if the field of view is greater than 180 • . In this case, Section 2.4 shows that the catadioptric projection functions p c and p + nc are not continuous at the infinite plane since two different limits are obtained. Thus, e i j is not a continuous function at the infinite plane (if t = 0). This error is not ideal since we do not have the C 2 continuity for the crossing of the infinite plane. Even worse, the C 0 discontinuity of e i j is very great since the two projection limits (if t converge to 0) are two antipodal points separated by the small circle of the catadioptric image. If X j goes from the right side (where ||e i j || is less than 2 pixels in practice) to the wrong side of the infinite plane, the resulting increase of ||e i j || may reach hundreds of pixels due to the great distance between antipodal projections (Section 2.3). Although useful for bundle adjustment convergence in our context, such plane crossings are difficult for an LM iteration which should decrease the cost function E (Section 3.1).

We see that infinite plane problems cannot be ignored for error e i j in our context: catadioptric cameras with fields of view greater than 180 • , scenes with distant 3D points, approximate calibration knowledge and image noise. In particular, point parametrization (x, y, z) is not sufficient and the use of the standard 2D reprojection error is not straightforward for bundle adjustment. A condition for an ideal error is C 2 continuity with 3D point parametrization (x, y, z, t), even at the infinite plane.

Image Error

This part presents an image error based on the standard 2D reprojection error of a catadioptric camera. We aim to improve the smoothness of this error for the crossing of the infinite plane (Section 3.2). Let X j be the four homogeneous coordinates of the (finite) j th point in the world coordinates. If the model is central, we define

e i+ j = p i c (X j ) -u i j , e i- j = p i c (-X j ) -u i j (8)
where p i c (.) and p i c (-.) are the antipodal projection functions of the i th camera (defined by Equations 1 and 2). If the model is non-central, we define

e i+ j = p i+ nc (X j ) -u i j , e i- j = p i- nc (X j ) -u i j ( 9
)
where p i+ nc and p i- nc are the antipodal projection functions of the i th camera (defined by Equations 4 and 5).

At first glance, the error function e i j might be defined by the standard 2D reprojection error: e i j = e i+ j . However, this function is discontinuous at the infinite plane (Section 3.2). We propose

e i j = e i+ j , if ||e i+ j || < ||e i- j || or if e i- j is not defined e i- j , if ||e i- j || < ||e i+ j ||
or if e i+ j is not defined [START_REF] Harris | A Combined Corner and Edge Detector[END_REF] and see that ||e i j || = min(||e i+ j ||, ||e i- j ||) if both e i+ j and e i- j are defined.

First, we show that ||e i+ j || = ||e i- j || is very improbable in practice during the minimization of Equation 6. This would imply values of ||e i j || = ||e i+ j || = ||e i- j || greater than the radius of the small circle due to the large distance between two antipodal projections (Section 2.3). However, all initial ||e i j || are small (less than 2 pixels) and the bundle adjustment attempts to decrease all ||e i j || simultaneously by minimizing their sum of squares.

Second, we show that e i j is C 2 continuous at finite point X j . If ||e i+ j || < ||e i- j || at X j , this inequality still holds in a neighborhood of X j thanks to the continuity of the antipodal functions. Thus, e i j = e i+ j in this neighborhood thanks to Equation 10. We see that e i j is C 2 continuous at the finite point X j just like p i+ nc and p i c . The proof is similar if

||e i- j || < ||e i+ j || at X j .
Third, we show that e i j is C 0 continuous at infinite point X ∞ . Let X j be a finite point which converges to X ∞ . Section 2.4 shows that the antipodal projections of X j converge to the projections u + ∞ and u - ∞ of X ∞ in the i th image. If both antipodal projections exist, e i+ j and e i- j converge to u + ∞u i j and u - ∞u i j . Thus, e i j converges to the vector among u + ∞u i j and u - ∞u i j which has the smallest modulus (Equation 10): e i j is continuous at X ∞ . If only one antipodal projection exists and converges to u ∞ , the corresponding error e i+ j (or e i- j ) converges to u ∞u i j and the other error e i- j (or e i+ j ) is not defined. Thus, e i j is defined by e i+ j (or e i- j ) and converges to u ∞u i j : e i j is continuous at X ∞ . If no antipodal projection exists, e i j is not defined at X ∞ .

Last, we show that e i j is the standard 2D reprojection error e i+ j if X j is finite and is on the right side of the i th mirror during minimization (i.e. if X j has not crossed the infinite plane). In this case, ||e i+ j || is small and the large distance between antipodal projections provides a large ||e i- j ||. Thus, e i j = e i+ j .

Angular Error

We also introduce an angular error for a catadioptric camera: error e i j is such that ||e i j || is the angle between two rays. The first one is a back-projected ray obtained by the camera calibration applied to u i j . The second one is a ray corresponding to the reprojection by the i th camera of the j th point. Let d i j and D i j be the directions of the first and second ray, respectively.

If the model is central, the two rays start from the camera center t i of the i th camera. Let R i be its rotation matrix, and X j the homogeneous coordinate of the j th point (t i , R i , X j are in the world coordinate system). The definitions of d i j and D i j are straightforward in the camera coordinate system:

d i j = C -1 (u i j ) ||C -1 (u i j )|| , D i j = R i [I 3 | -t i ]X j ||R i [I 3 | -t i ]X j || . (11)
If the model is non-central, s i j , n i j and a i j are introduced. Point s i j is the intersection of the mirror surface and the back-projected ray by the perspective camera of u i j . We also define n i j as the mirror normal at s i j and a i j as the direction of the back-projected ray (pointing outside the mirror). In this noncentral context, we redefine R i as the orientation and t i as the origin of the i th mirror coordinate system in the world coordinate system. The first ray is a half line starting from s i j with the direction d i j defined by the reflection law. Both s i j and d i j are expressed in the i th mirror coordinates, and are fixed by u i j and all parameters of the perspective camera. The second ray is the half line starting from s i j towards X j . The expressions of d i j and D i j are easy given s i j , n i j , a i j and X t j the fourth homogeneous coordinate of the j th point:

d i j = 2(n i j a i j )n i j -a i j ||2(n i j a i j )n i j -a i j || , D i j = R i [I 3 | -t i ]X j -X t j s i j ||R i [I 3 | -t i ]X j -X t j s i j || . ( 12
)
Now, the main subject in this part is the choice of function e i j with the C 2 continuity (Section 3.1) even if X j is at the infinite plane (Section 3.2). At first glance, we might choose the 1D error e i j = arccos(d i j .D i j ). Unfortunately, Appendix E shows that this function is not C 1 continuous if e i j = 0. This C 1 discontinuity at the exact solution complicates the convergence in practice. A second try might be e i j = f (d i j .D i j ) with f a decreasing C 2 continuous function such that f (1) = 0, if we accept that e i j is not an angle. The resulting e i j is C 2 continuous and e i j ≥ 0. We deduce that e i j has a local extrema if e i j = 0: the Jacobian of e i j is zero here. The convergence rate of LM is reduced in this context (Section 3.1).

Our final e i j proposition does not have these problems and is defined by

e i j = π 2 (R i j D i j ) with R i j a rotation such that R i j d i j =    0 0 1    . ( 13
)
Projection π 2 was defined in Equations 3. Since ||π 2 (x, y, z

)|| = x 2 +y 2 z 2
is the tangent of the angle between ( 0 0 1 ) and ( x y z ) , we see that

||e i j || 2 = tan 2 (α i j ) with α i j = angle(d i j , D i j ).
This result is independent of the choice of R i j . The resulting cost function E (defined in Equation 6) is a sum of squared tangents of angles between rays. Error e i j is well defined (i.e. angle(d i j ,

D i j ) = π 2 [π]
) and E is a good approximation of the sum of the squared angles in our context since these angles are small in practice. Furthermore,

e i j = π 2 (R i j R i [I 3 | -t i -R i s i j ]X j ) (14)
since the projection π 2 cancels the scale of D i j defined in Equations 11 and 12 (s i j = 0 in the central case). We recognize a projection by a perspective camera such that the orientation and center are parametrized. So, e i j inherits all high smoothness properties of the perspective camera: e i j is C 2 continuous even if X j is in the infinite plane. This smoothness is obvious if the calibration is not refined by the BA. If it is refined by BA, we should assume that the mirror profile is such that d i j (and s i j in the non-central case) are C 2 continuous.

Overview of the Automatic Method

Now, the automatic Structure from Motion method is described for both catadioptric camera models using bundle adjustments (BA) proposed in Sections 3.3 and 3.4. First, the geometry of the image sequence is estimated for the central model. Then, a second estimation using the non-central model is obtained. In both cases, an approximate and given calibration is refined with the 3D. Many technical details (useful for possible re-implementers) about matching and initialization steps are also given in the Appendix.

Assumptions for Calibration Initializations

The following assumption is required: a surface-of-revolution mirror, whose lower and upper circular cross sections are visible, is placed in front of a natural perspective camera (zero skew, aspect ratio set to 1). Furthermore, the perspective camera must point towards the mirror and the perspective center is in the immediate neighborhood of the mirror symmetry axis.

In this context, the projections of the two mirror circular cross sections are approximated by two circles in images. The central model only requires the knowledge of these two circles and the radially symmetric approximation.

The non-central model only requires the knowledge of the mirror profile. The assumptions above are needed for the calibration initialization steps in Sections 4.2 and 4.5. These steps are not the main subjects of this paper, and more general methods are available if we assume that the projections of the two mirror circular cross sections are general ellipses (e.g. [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF] in the central case and [START_REF] Wu | Camera Calibration from the Quasi-affine Invariance of Two Parallel Circles[END_REF] in the non-central case).

Central Calibration Initialization

First, the large and small circles in each catadioptric image are detected and estimated using RANSAC and Levenberg-Marquardt methods applied on the vertices of regularly polygonized contours. Then, the central calibration C (Equation 2) is initialized as follows. We define r(α) by the linear function such that r(α up ) = r i up and r(α down ) = r i down . The image circle radii r i down , r i up are obtained in the circle estimation step, and the field of view angles α down , α up are given by the mirror manufacturer. These angles are not exactly known since they depend on the relative position between the mirror and the pinhole camera.

Central Reconstruction Initialization

Harris points are detected and matched for each pair of consecutive images in the sequence using correlation without any epipolar constraint. The corresponding ray directions are also obtained from the central calibration initialization. Then, the essential matrices for these pairs are estimated by RANSAC (using the 7 point algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]) and refined by Levenberg-Marquardt. 3D points are also reconstructed for each pair. Many of these points are tracked in three images, and they are used to initialize the relative 3D scale between two consecutive image pairs. Now, points in 3 views are reconstructed and we obtain the reconstruction of each triple of consecutive images. More details about these matching and estimation steps are given in Appendix A and B.

Last, the full sequence geometry is obtained by many BAs applied in a hierarchical framework to merge all partial geometries [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Once the geometries of the two camera sub-sequences 1

• • • n 2 , n 2 + 1 and n 2 , n 2 + 1, • • • n are
estimated, the latter is mapped in the coordinate system of the former thanks to the two common cameras n 2 , n 2 + 1, and the resulting sequence 1 • • • n is refined by a BA. The angular error (Equations 11 and 13) is preferred for these BAs since we found its convergence more robust in practice. A final BA with image error (Equations 8 and 10) is applied to the full sequence to obtain a MLE with the common Gaussian model in images (Section 1.3).

Central Reconstruction and Calibration Refinements

Once the full central geometry is obtained for the approximate (linear) function r(α) defined above, r(α) is redefined as a cubic polynomial whose the 4 coefficients should be estimated. We assume that the central calibration is constant and apply an additional BA using image error (Equations 8 and 10) to estimate the 4 + 6c + 3p parameters of the sequence (c is the number of cameras, p is the number of 3D points).

Non-Central Calibration Initialization

The parameters K p , R p and t p of the perspective camera (Section 2.2) are initialized from the detected large and small circles. These circles are the images of the known mirror circular cross sections with radii r up and r down in z-plane z = z up and z = z down , respectively. From Section 4.1, we have R p ≈ I 3 , t p = ( -x p -y p -z p ) where max(|x p |, |y p |) z p , and K p is the 3 × 3 diagonal matrix with focal length f p .

First, an approximate value of z p is obtained by measuring the distance between the mirror and the perspective camera. Second, f p is estimated assuming x p = y p = 0 by the Thales relation r i up /f p = r up /(z p +z up ). Third, R p and x p , y p are estimated by projecting the circular cross section centers. More details are given in Appendix C.

Non-Central Reconstruction Initialization

The reconstruction with the non-central model is initialized from the reconstruction with the central model. A first possibility is the approximation of each non-central camera by a central camera such that the center is the mirror apex: the i th mirror coordinate system is defined by the pose (t i , R i ) of the i th central camera, and the j th 3D point X j is retained as such. The 3D scale factor of the scene, in fact, is not a free parameter as in the central case since the non-central image projections are dependent on it. We choose the initial 3D scale factor by multiplying t i and X j by λ ∈ such that

d = 1 c-1 c-1 i=1 ||λt i -λt i-1
|| is physically plausible. The value of d is obtained by an estimate of the step length between two consecutive images.

Then, a BA is applied to refine the 6c + 3p parameters of the non-central reconstruction. We found that the angular error (Equations 12 and 13) is better than the image error (Equations 9 and 10) to start parameter refinement: the convergence is more robust and one LM iteration is about 3 times faster. Last, we finish parameter refinement by applying a BA with the image error to obtain a MLE of the geometry with the common Gaussian model in images (Section 1.3).

Non-Central Reconstruction and Calibration Refinements

Once the full non-central geometry is obtained from the method above, additional parameters among K p , R p and t p are selected to be refined. We assume that t p and f p are constant during the camera motion in the scene, but the orientation R p is perturbed around I 3 (slight rotations of the perspective camera are possible around a screw). Thanks to the mirror symmetry around the z-axis (Section 2.2), R p has only two parameters θ x and θ y : R p = R x (θ x )R y (θ y ) with rotations R x and R y around x and y-axis of the mirror coordinate system. The number of additional parameters is k = 1 + 3 + 2c. Finally, an additional BA with image error (Equations 9 and 10) is applied to refine the k + 6c + 3p parameters of the sequence.

Experiments

After the description of the experimental context in Section 5.1, specific results are given for the central model in Section 5.2 and for the non-central model in Section 5.3. Both models are compared with many criteria including accuracy and uncertainty in Section 5.4. Section 5.5 compares reconstructions obtained by minimizing angular and image errors. A piecewise planar 3D model and some technical details are also given in Sections 5.6 and 5.7.

Experimental Context

The user moves along a trajectory on the ground with the omnidirectional system mounted on a monopod, alternating a step forward and a shot. We use a cheap system designed for panoramic picture generation and imagebased rendering from a single view point, given a single shot of the scene (see Figure 5). The mirror is not a quadric and has a known profile (Equation 17) such that we have an approximately uniform resolution in the image radial direction (r(α) of Equation 2 is linear). Such a catadioptric camera/system is called an "equiangular" camera, it is not a central camera, and the size of the caustic profile [START_REF] Grossberg | A general imaging model and a method for finding its parameters[END_REF] is about half the size of the mirror profile (right of Figure 5). The symmetry axes of camera and mirror are not exactly the same, since the axes alignment is manually adjusted and visually checked.

Central Model

Panoramic images obtained from one image of the tested sequences are shown in Figure 6. Top views of resulting reconstructions by methods described in Sections 4.2, and 4.3 are shown in Figure 7 and 11. These results are obtained with a linear calibration function r(α) defined by the field of view angles α up = 40 • , α down = 140 • given by the mirror manufacturer.

Fountain

The Fountain sequence is composed of 38 images of a background city and a close-up fountain at the center of a traffic circle. We have found that about 50% of recovered essential matrices are obviously incorrect for this sequence, since the epipoles are roughly orthogonal to the camera motion as shown in Figure 8. However, the 3-view calculations remove all blunders of 2-view calculations thanks to the 3-view selection of 3D points (Appendix B) and angular bundle adjustment with inlier update (Section 5.7). The recovered camera motion is smooth and circular around the fountain. This result is consistent with our knowledge of the true camera motion.

Both unclosed and closed versions of this sequence are reconstructed. The unclosed version is obtained by duplication at the sequence end of the first image, and is useful since it provides information about the pose accuracy by measuring the gap between both sequence ends (this gap is significant of the drift of the unclosed reconstruction process). The distance ||t 0t 38 || between both sequence ends is 0.01 times the trajectory diameter, and the angle of the relative orientation R 38 (R 0 )

-1 is 0.63 • . The closed version is considered in the rest of the paper (including Figures and Tables).

Road and House

The Road sequence is composed of 54 images taken along a little road on flat ground. The background includes buildings, parking and fir trees. All recovered essential matrices seem to be correct according to the positions of the epipoles. The House sequence is composed of 112 images taken in a cosy house, starting in the living room, crossing the lobby, having a loop in the kitchen, re-crossing the lobby and entering a bedroom.

Robustness to the given Field of View Angles

This section shows the robustness of the central reconstruction methods to linear calibration function r(α) defined by rough values of α up and α down (Section 4.2). This is useful in practice, since the field of view angles are sometimes unknown or inaccurate (they depend on the relative position between the mirror and the pinhole camera). We also experiment the r(α) refinement.

The experiments are summarized in Figure 9. For each r(α) initialization defined by (α up , α down ) ∈ {(40, 140), [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF]160), (60, 120), [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF]120), (60, 120)}, the central methods are applied. The reconstruction initialization (Section 4.3) fails twice for the Fountain with these large inaccuracies of ±20 • . No failure occurs for the easier case ±10 • . The last step is the simultaneous reconstruction and calibration refinements (Section 4.4). Function r(α) is redefined as a cubic polynomial and is estimated using bundle adjustment: the four polynomial coefficients are new unknowns, and the set of inliers is updated four times during the minimization. The recovered r(α) are shown on the right of Figure 9 for each initial (α up , α down ) and each of the 3 sequences. The exact r(α) does not exist since the catadioptric camera is non-central: it depends on the depth of points. It is also dependent on the settings of the perspective camera which are slightly different for the 3 sequences. However, we assume that the expected results are near the linear r(α) defined by the manufacturer angles. With this in mind, we see that the calibration improvement is significant. Furthermore, the r(α) refinements are stable since the calibration curves are very similar for each sequence, except in the neighborhood of the small circle (radius r i down ) for the Fountain images. This default is probably due to the high lack of image matches that we observed near the small circle in the whole sequence. The means of recovered (α up , α down ) are (46 • , 149 • ), (44 • , 139 • ), (45 • , 144 • ) for Fountain, Road and House. These angles are slightly greater than the given manufacturer angles (40 • , 140 • ).

Non-Central Model

Non-central methods are also applied to enforce the mirror knowledge in the geometry estimation (the mirror profile is defined by Equation 17). Calibration and reconstruction are initialized thanks to estimates of z p and the 3D scale factor of the scene as mentioned in Section 4.5 and Section 4.6, respectively. The largest z p available with the monopod is chosen to increase the depth of field and obtain the entire mirror in sharp focus with the perspective camera. We measure z p = 48 cm between the mirror and the camera. The approximate trajectory lengths of Fountain, Road, and House sequences are 16, 42 and 22 meters, respectively. The non-central reconstructions are qualitatively similar to the central ones in Figures 7 and11. The experiments show that the refinements of calibration and 3D scale are difficult in our context. The reasons are given in Section 5.3.1 for the calibration and in Section 5.3.2 for the 3D scale factor.

Calibration Refinements

As mentioned in Section 4.7, there are 1 + 3 + 2c parameters of the perspective camera to be refined: one focal length f p , one center t p = ( -x p -y p -z p ) and many rotations R p = R x (θ x )R y (θ y ). Obviously, these parameters are added to the 6c + 3p parameters of the 3D in the bundle adjustment.

We note that the mirror size (radius 3.7 cm and axis length 3.5 cm) is small in comparison with the distance between the pinhole center and mirror (z p = 48 cm). In this context, small perturbations of θ x , θ y and f p are almost compensated by certain perturbations of t p to remain the non-central unchanged projection of a 3D point. More precisely, we have a rotationtranslation ambiguity [START_REF] Adiv | Inherent ambiguities in recovering 3-D motion and structure from a noisy flow field[END_REF]: it is difficult to refine R p and the x-y components of t p simultaneously. Also, we have ambiguity between zoom and translation small scale scene (r t = 25 cm) large scale scene (r t = 2.5 m) 1 The ground truth reconstructions are perturbed by initial homothety s init and image Gaussian noise of σ = 1 pixel. For each s init , a homothety s ref ined is estimated by non-central bundle adjustment. Homothety s ref ined is the ratio between the recovered scale factor and its true value (the ideal result is s ref ined = 1). The small scale scene has the best values of s ref ined .

s init 0.5 .
along the focal axis [START_REF] Bougnoux | From Projective to Euclidean Space under any practical situation, a criticism of self-calibration[END_REF]: it is difficult to simultaneously refine f p and the component of t p along the third column of R p . Obviously, these ambiguities are inherent to the estimation problem (a different method will not remove them).

In the experiments, convergence is slow and the final f p , t p , R p values are similar to their initial values. The ambiguities are confirmed by high correlation coefficients for our sequences (obtained from the covariance matrix of the estimated parameters [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]): we always have |σ fp,zp | ≥ 0.999, and |σ θx,yp | ≥ 0.8, |σ θy,xp | ≥ 0.8 in the majority of cases.

3D Scale factor Estimation

Refinement of the 3D scale factor of the scene by bundle adjustment is theoretically possible using a non-central catadioptric system, since the non-central image projection changes with the scale factor. However, experiments on our image sequences show that the scales recovered by the non-central reconstruction initialization and refinement (Sections 4.6 and 4.7) are similar to their initial values. The reason is the following: a scale change of the scene (including the mirror trajectory) does not lead to significant changes to the projections of any 3D point distant enough from the mirrors.

The following synthetic experiment shows that it is more difficult to estimate the 3D scale factor for large scale scenes than small scale scenes. Two ground truth reconstructions are defined by half turn of a Fountain-like scene with trajectory radius r t , mirror orientations R i perturbed around I 3 , 20 camera poses and 1000 points well distributed in 3D and 2D spaces. The only 3D difference between both reconstructions is the 3D scale of the scene points and the mirror trajectory defined by r t = 2.5 m and r t = 25 cm. The perspective parameters K p , R p , t p are the same and the (exact) image projections are different. First, image projections are corrupted by a Gaussian noise of σ = 1 pixel, and all mirror locations t i and points X j are multiplied by an initial factor s init . Second, the bundle adjustments of Section 4. [START_REF] Faugeras | The Geometry of Multiple Images[END_REF] 2 The RMS in pixels, the numbers of 3D reconstructed and 2D detected points (inliers) for four reconstructions: central reconstructions with initial and refined calibrations, non-central reconstructions with initial and refined calibrations. The degree of freedom (d.o.f.) depends on the numbers c and p of cameras and 3D points.

to these perturbed reconstructions enforcing the exact perspective parameters (z p = 48 cm) and taking into account the outliers. Table 1 shows for many values of s init the resulting RMS (pixels) and ratio s ref ined between the recovered scale factor and its exact value. We note that the RMS has no clear minimum for the large scale scene with r t = 2.5 m, and that the scale factor estimations are best for the small scale scene with r t = 25 cm.

In practice, we have abandoned attemps to estimate the 3D scale factor accurately: the measure of distance between two scene points or camera centers is too difficult with this catadioptric system and usual scenes.

Performance Comparisons Between Central and Non-Central Models

Quantitative comparisons between real 3D reconstructions using the central and non-central models are given.

Consistency

Table 2 shows consistencies between many reconstructions and the multi-view matching for the sequences Fountain, Road and House. There are four reconstructions: central reconstructions with initial (Section 4.3) and refined (Section 4.4) calibrations, non-central reconstructions with initial (Section 4.6) and refined (Section 4.7) calibrations. The consistency criteria are the RMS, the numbers of 3D and 2D points which are consistent (3D reconstructed and 2D detected points X j and u i j are consistent if the standard 2D reprojection error ||e i j || is less than 2 pixels). The non-central reconstructions have the best consistencies: improvements about 5% (sometimes 12% for the Road) are obtained for the numbers of 3D and 2D points, with slightly slower RMS. We also see that the consistency improvements by calibration refinement are non negligible for the central model. They are negligible for the non-central model.

Difference

Two reconstructions "a" and "b" are compared as follows. The camera location difference and 3D point difference are respectively

E (a,b) t = 1 I i ||S(t i b ) -t i a || 2 , E (a,b) x = 1 J j ||S(X b j ) -X a j || 2 ||t i(j) a -X a j || 2 (15)
with S the similarity transformation minimizing E , I the number of cameras, J the number of 3D points, and i(j) the index of the closest t i a to the j th point X a j . t i a is the apex of the i th mirror in the non-central case and the i th camera center in the central case.

The 3D differences between central (t i c , X c j ) and non-central (t i nc , X nc j ) reconstructions are the followings. We obtain E

Pose accuracy

A real sequence (Figure 10) is taken in an indoor controlled environment: the motion of the catadioptric system is measured on a rail, in a 7m × 5m × 3m room. The trajectory is a 1 meter long straight line by translation, with 6 equidistant and aligned poses. The location errors are E (g,c) t and E (g,nc) t with t i g the location ground truth and E (.,.) t defined in Equations 15. We obtain E (g,c) t = 1.1 mm and E (g,nc) t = 1.2 mm for the central and non-central models, respectively. Both models provides similar and good accuracies for the location estimation in this context. The angle E r = max i,j {arccos( 1 2 (trace(R i (R j ) -1 ) -1))} is our orientation error since the ground truth of R i is unknown and maintained constant. The results are E c r = 0.5 3 Right: the camera center and 3D point uncertainties (cm) for central reconstructions. Left: the lengths (m) of the trajectories and the gauge constraints. The uncertainties are the length of the major semi-axis of the uncertainty ellipsoids for the probability 90%. u 3. Ellipsoids are very long for 3D points which are distant or roughly aligned with the cameras from which they are reconstructed.

Uncertainty

When ground truth is not available, informations about the reconstruction quality are provided by the uncertainty (covariance matrix) estimation of the geometry parameters [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. This estimation requires the common Gaussian model (Section 1.3) for the image error minimized by bundle adjustment.

A trivial camera-based gauge constraints is chosen to obtain minimal parameterizations [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. Since the central reconstruction is defined up to a similarity transformation (7 d.o.f.), we choose the constraints R 0 = I 3 , t 0 = 0 and one fixed coordinate of t i 0 . The resulting uncertainties for 3D points and camera centers are given in Table 3 (also refer to Figure 11 for the Road sequence). In the non-central case, only R 0 = I 3 , t 0 = 0 should be sufficient since the global scale of the scene is fixed by the mirror size. In practice, we found that the medians of the resulting non-central uncertainties are 3-26 greater than the medians of the central uncertainties. The reason is given in Section 5.3.2: the 3D scale factor of the scene cannot be estimated accurately.

Angular vs. Image Error

The minimization of angle errors (Section 3.4) may appear to be a heuristic choice. It is shown in this part that the results are not so different to those obtained with the more standard minimization of image errors (Section 3.3). We use the ground truth reconstructions introduced in Section 5.3.2 and corrupt the image projections by a Gaussian noise of σ = 1 pixel. Then, we compare the reconstructions (t i A , X A j ) and (t i I , X I j ) obtained by a BA minimizing angular and image errors, respectively with the ground truth (t i g , X g j ). The BAs are initialized by the noisy ground truth reconstruction, the calibration and inliers are maintained to be the same, and the reconstruction comparisons are done with the measures defined by Equations 15.

The results are E (g,A) t = 0.29 cm, E (g,A) x = 0.027, E (g,I) t = 0.28 cm and E (g,I) x = 0.027 for the non-central and large scale reconstruction of Section 5.3.2. They are similar with the central model and/or the small scale reconstruction.

Piecewise Planar 3D Modeling

For each planar piece of the targeted 3D model, we choose a catadioptric image of the sequence and define the contour of the piece manually. Then, we select reconstructed points by their projections in the delimited region and estimate the support plane from these points by minimizing a sum of (squared) pointto-plane distances. The Mahalanobis point-to-plane distance [START_REF] Schindler | On Robust Regression in Photogrammetric Point Clouds[END_REF] is preferred to the Euclidean distance to favor the points with the least uncertainties. In this context, independent point covariances are useful in order to obtain a Maximum Likelihood Estimation of the plane. So these covariances are estimated by inverting each Hessian of the independent ray intersection problems.

A piecewise planar model of the House is shown in Figure 12. It is not difficult to guess where the living-room, the kitchen, the lobby, and the bedroom are. 36 planes are estimated from the non-central reconstruction. Plane accuracy depends on the number of 3D points selected (between 7 and 178), their accuracies, uncertainties, and distributions in images.

Some Technical Details

Computation times for 1632 × 1224-images with a P4 2.8GHz/800MHz are about 10 seconds for each single image calculation (points and edge detections, circles), 20 seconds for each image pair calculation (matching by quasi-dense propagation [START_REF] Lhuillier | Match Propagation for Image-Based Modeling and Rendering[END_REF], essential matrix estimation) for each sequence, and 5 minutes for the hierarchical reconstruction step applied to the House sequence (total time: 1 hour). About 700-1100 point matches satisfy the epipolar constraint between two consecutive images.

An image point is considered as an outlier for all angular (respectively, image) bundle adjustments if the corresponding error is greater than 0.04 radians (respectively, 2 pixels). The sets of inliers and outliers are updated one time during each bundle adjustment. Accuracy increases when new inliers are discovered and involved in the score to minimize. Furthermore, inliers which becomes outliers after the first round of bundle adjustment are ignored in the second round to gain robustness. Usual and final RMS are about 0.005-0.006 radians (respectively, 0.7-0.8 pixels) with RMS= E/(#2D). Function E is defined by Equation 6 and #2D is the number of e i j involved in E.

All bundle adjustments are implemented using fully analytical derivatives, except the image error in the non-central case which requires an iterative method for M + and M -calculations (more details in Appendix D). One bundle adjustment iteration for the final refinement of the House takes about 1.6 seconds using the non-central image error (7 seconds are also required for M + and M -initialization before all optimizations), and about 0.6 seconds for others bundle adjustments.

Conclusion

The methods described in this paper provide an automatic, robust and optimal estimation of the scene structure and camera motion for image sequences acquired by a catadioptric camera. First, we propose bundle adjustments minimizing angle and image errors, by taking care of the targeted smoothness conditions for good convergence. The image error provides a Maximum Likelihood Estimate of the sequence geometry for the common Gaussian model in images, and its smoothness is improved in the infinite plane. The angle error has the ideal smoothness (C 2 continuity, even at the infinite plane). The second contribution is an extensive experimental study in a context that we expect representative and useful for applications: a hand-held and equiangular catadioptric camera moving on the ground. Many experiments are presented about robustness (including the initialization robustness to the given field of view angles), accuracy and uncertainty estimations, performance comparisons between central and non-central models, and piecewise planar 3D modeling from the reconstructed points. The central model provides good approximation of the (real) non-central model. On the other hand, the 3D scale factor is difficult to estimate with the non-central model. The only demonstrated improvement by the non-central model is obtained for the consistency between matching and geometry. We also discuss calibration refinements for both central and non-central models. Last, the 3D reconstruction system is described as a whole.

Future works include applications (image-based modeling and rendering, vehicle localization ...), matching improvements, and camera calibration by enforcing constraints on the 3D scene (especially for the non-central case which has rarely been discussed until now).

portion of the distortions is compensated for by image rotation around the circle center. The Harris point detector [START_REF] Harris | A Combined Corner and Edge Detector[END_REF] is used because it is invariant to such rotations and it has good detection stability. We also compensate for the rotation in the neighborhood of the detected points before comparing the luminance neighborhood of two points using the ZNCC score (Zero Mean Normalized Cross Correlation). To avoid incorrect matching due to repetitive textures, the following procedure is applied. First, the points of interest of an image are matched to other points of interest in the same image. The result is a "reduced" list of points which are not similar to others according to ZNCC in their corresponding search area. Second, the points of the reduced lists of two different images are matched applying the same correlation score, search areas and thresholds. Now the matching errors due to repetitive textures have been greatly reduced, but the current list of matches is very incomplete. Third, this list has been completed thanks to a quasi-dense match propagation [START_REF] Lhuillier | Match Propagation for Image-Based Modeling and Rendering[END_REF]: the majority of image pixels are progressively matched using a 2D-disparity gradient limit and the uniqueness constraint, and two interest points roughly satisfying the resulting correspondence mapping between both images are added to the list. The window sizes and lower bound thresholds for ZNCC are the same as in [START_REF] Lhuillier | Match Propagation for Image-Based Modeling and Rendering[END_REF]. The resulting list of matched interest points is used in the epipolar geometry estimation step (described in Appendix B).

Appendix B. Two-View and Three-View Central Initialization

More details are given on the geometry initialization described in Section 4.3. Let d 0 j and d 1 j be the ray directions (normalized vectors) obtained by the central calibration applied to the j th pair of matched points in images 0 and 1. First, a fundamental matrix F is obtained with the 7-point algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] : (1) each pair (d 0 j , d 1 j ) provides a normalized linear equation for the 9 parameters of F, (2) F = F 1 + λF 2 where F 1 , F 2 are two solutions of a 7 × 9 linear system and (3) F is obtained by solving a cubic polynomial equation in λ to enforce the constraint det(F) = 0. Second, an essential matrix E is obtained from F by forcing the two largest singular values of F to be the same by SVD [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Third, we count the number of pairs (d 0 j , d 1 j ) such that the angle between d 1 j and the normal Ed 0 j ||Ed 0 j || of the epipolar plane of d 0 j is equal to π 2 up to a threshold t 0 . This process is repeated many times by the RANSAC method, and the E with the largest number of pairs is retained. The E refinement by Levenberg-Marquardt is straightforward: we use the parametrization E(t, R) = [t] × R and minimize a Longuet-Higgins criterion [START_REF] Faugeras | Three-Dimensional Computer Vision -A Geometric Viewpoint[END_REF] defined by LH(t, R) = j (d with x(α), y(α), z(α) three real C 1 continuous functions with parameter α such that ( x(0) y(0) z(0) ) = ( 0 0 1 ). Now, we show that the limit of ∂e i j ∂α is not well defined if α converges to 0.

Fig. 2 .

 2 Fig.2. Left: two concentric circles of radii r i up , r i down in the image. Right: angles α up , α down , α from the z-axis define the field of view for a point X by α up ≤ α(X) ≤ α down . The circle "up" (respectively, "down") of rays on the right is projected by the central model on the large (respectively, small) circle on the left.

Fig. 4 .

 4 Fig. 4. This figure is obtained if the finite point X of Figure 3 converges toward the right to the infinite point X + ∞ . The two lines which go across M + (t p , X + ∞ ) and M -(t p , X + ∞ ) have the same direction X + ∞ .

Fig. 5 .

 5 Fig. 5. Left: the 360 One VR (Kaidan) mirror with the Coolpix 8700 (Nikon) camera, mounted on a monopod. Middle: the view field is 360 • in the horizontal plane and about 50 • above and below if the camera is pointing toward the sky. Right: the profile of the mirror caustic for z p = 48 cm and x p = y p = 0.

Fig. 6 .

 6 Fig. 6. Panoramic images obtained from one omnidirectional image of sequences Fountain, House, and Road. They are given to help understand the scenes, but they are not used by the methods.

Fig. 7 .

 7 Fig. 7. Top views of Fountain (38 views, 5857 points) and House (112 views, 15504 points) central reconstructions. The Road (54 views, 10178 points) reconstruction is shown in Figure 11. Several of these points are difficult to reconstruct accurately if they are distant or roughly aligned with the cameras from which they are reconstructed. These results are similar for central (with linear and cubic r(α) functions) and non-central models.

Fig. 8 .

 8 Fig. 8. Epipolar curves for images of two estimated pair-wise central geometries of the Fountain sequence, and true camera motion indicated by the black arrows. The geometry estimation is obviously incorrect on the right since the epipoles do not agree with the true camera motion (they do on the left). Such blunders are corrected by the three view calculations.

Fig. 9 .

 9 Fig. 9. From left to right: many initial α(r), and the resulting refined α(r) for the Fountain, Road and House sequences, respectively. Function α(r) is the reciprocal function of the radial function r(α) involved in the central calibration (Equation2). The ranges are always [r i down , r i up ] × [0, π]. The r-distribution of inlier matches {m} is also given by histogram for each sequence, at the bottom.

Fig. 10 .

 10 Fig. 10. A panoramic image from the "controlled" sequence.

(

  nc,c) t = 0.79 cm, E (nc,c) x = 0.027 for the Fountain (respectively, E (nc,c) t = 2.6 cm, E (nc,c) x = 0.017 and E (nc,c) t = 2.75 cm, E (nc,c) x = 0.054 for the Road and House).

  are respectively the rank 0/4 (smallest), rank 1/4, rank 2/4 (median), rank 3/4 and rank 4/4 (largest) semi-axis lengths for the 3D points. u c is the largest semi-axis length for the cameras.

Fig. 11 .

 11 Fig. 11. Top view of the Road central reconstruction (54 views, 10178 points) and the uncertainty ellipsoids of Table3. Ellipsoids are very long for 3D points which are distant or roughly aligned with the cameras from which they are reconstructed.

Fig. 12 .

 12 Fig. 12. A piecewise planar model of the House.

j

  || )2 . The last step of the two-view geometry initialization is the reconstruction of each point X j by minimizing the function X j → ||e 0 j || 2 + ||e 1 j || 2 with e i j defined by Equations 11 and 13.

Fig. 13 .

 13 Fig. 13. Left: useful notations for non-central calibration initialization assuming x p = y p = 0. Right: t p is located in a circle in the plane [z = -z p ] given z p and the angle between directions from t p toward the centers of large and small border circles.

  are applied

	Calib.	initial central	initial non-central	refined central	refined non-central
	Criteria	#3D,#2D,RMS	#3D,#2D,RMS	#3D,#2D,RMS	#3D,#2D,RMS
	Fountain	5857,31028,0.84	6221,33262,0.78	5953,32214,0.77	6223,33876,0.76
	Road	10178,50225,0.87 11312,56862,0.82	10814,55297,0.77	11313,57148,0.79
	House	15504,75100,0.83 16432,80169,0.78	15800,77883,0.76	16447,80311,0.76
	d.o.f. Table	6c + 3p -7	6c + 3p -6	4 + 6c + 3p -7	4 + 8c + 3p -6

  • and E nc r = 0.35 • .

		Length Gauge constraints	u c	u	0/4 p	u	1/4 p	u	2/4 p	u	3/4 p	u	4/4 p
	Fountain	16 m	fixed R 0 , t 0 , t 14 x	1.02	22 50.7 102 230 27e+4
	Road	42 m	fixed R 0 , t 0 , t 53 x	4.60 4.37 7.00 23.1 130 44e+4
	House	22 m	fixed R 0 , t 0 , t 111 x	7.62 5.17 7.27 12.2 26.7 11e+3
	Table

Appendix A. Two-View Matching

The usual matching procedure of interest points using correlation cannot be directly applied to the omnidirectional images for two reasons: (1) the matching ambiguity due to the repetitive textures in one image and [START_REF] Aliaga | Accurate Catadioptric Calibration for Real-time pose estimation in Room-size Environments[END_REF] the geometric distortions between matched patches of two images taken from different view points. Previously published matching methods deal only with one of these problems (e.g. relaxation [START_REF] Zhang | A Robust Technique for Matching Two Uncalibrated Images through the Recovery of the Unknown Epipolar Geometry[END_REF] is used for repetitive textures, regions [START_REF] Matas | Robust wide baseline stereo from maximally stable extremal regions[END_REF] or appropriate correlation windows [START_REF] Svoboda | Matching in Catadioptric Images with Appropriate Windows and Outliers Removal[END_REF] for geometric distortions in omnidirectional images).

In the context of the usual camera motions (roughly, translation motions with the pinhole camera pointing toward the sky), we observe that a high pro-Once the geometries of image pairs (0, 1) and (1, 2) are estimated with the method above, we know the central poses (t 0 , R 0 ), (t 1 , R 1 ) and (λt 2 , R 2 ) of the geometry of the image triple (0, 1, 2) up to the relative 3D scale λ ∈ with t 1 = 0. The 1 point RANSAC below is used to estimate λ. For each X j detected in image triple (0, 1, 2) and reconstructed from image pair (0, 1), λ j is estimated by minimizing the angular error in image 2 defined by

This error is derived from Equation 14 with s i j = 0. Then, we choose the λ j with the greatest number of points X j such that E 2 j (λ j ) is less than a threshold t 1 .

The point X j is retained in the geometry of the image triple (0, 1, 2) if it is detected in images 0, 1 and 2 (2 views are not enough for robustness). Furthermore, X j is reconstructed by minimizing the angular cost function

and it should satisfy max(||e 0 j ||, ||e 1 j ||, ||e 2 j ||) ≤ t 2 with t 2 a threshold. Last, a bundle adjustment with angular error (Equations 11 and 13) is used to refine the complete geometry of the image triple. Only one angular threshold t 0 = t 1 = t 2 = 0.04 (radians) is used in practice.

Appendix C. Non-Central Calibration Initialization

More details on the calibration initialization described in Section 4.5 are given in this Appendix. We estimate K p = diag(f p , f p , 1), t p = ( -x p -y p -z p ) and R p in the mirror coordinate system. First, an approximate value of z p is obtained by measuring the distance between the mirror apex and the perspective camera. Second, f p is estimated assuming x p = y p = 0 by the Thales relation r i up /f p = r up /(z p + z up ) (left of Figure 13).

Let c i up and c i down be the projections by the perspective camera of the centers of the two mirror circular cross sections. Thanks to the hypotheses (Section 4.1), c i up and c i down are approximated by the known centers of the detected large and small circles. Since the angle between directions pointing toward the centers of mirror border circles is that between vectors K -1 p c i up and K -1 p c i down , we know the radius r t of the circle in the plane [z = -z p ] where t p is located (right of Figure 13). Any t p in this circle is possible by the over-parametrization (R p , t p , R) of the non-central model (Section 2.2). Now, R p and x p , y p are estimated by projecting the circular cross section centers:

These equations and the hypothesis R p ≈ I 3 imply

x p y p .

Since any ( x p y p ) is possible in the circle x 2 p + y 2 p = r 2 t , we can choose (

Appendix D. Estimation and Derivation of M + (A, B)

The bundle adjustment with the non-central image error (Equations 9 and 10) requires efficient estimation and differentiation of M + (A, B). This function gives the reflection point on the mirror surface for the ray which goes across points A and B. Once these computations are done at (t p , X) with pinhole center t p and scene point X (in the mirror coordinate system), all derivatives of the projection p + nc (X) (Equation 4) are analytical according to the Chain Rule. Calculations are very similar for M -(A, B) and p - nc (X) (Equation 5).

The mirror surface is defined by the cylindric parameterization

with the mirror profile z(r) = 0.0287r + 0.218r 2 -0.0156r 3 + 0.00537r 4 , 0 ≤ r ≤ 3.7 cm. [START_REF] Lu | Fast and Globally Convergent Pose Estimation from Video Images[END_REF] Given A, B ∈ 3 , (r, θ) is estimated such that f (r, θ) is the reflection point M + (A, B) onto the mirror surface. Thus we have

by the reflection law, with n(r, θ) the mirror normal at point f (r, θ). Only 2 among these 3 equations are independent: (r, θ) is estimated such that g(r, θ, A, B) = 0 where

and R z (-θ) is the rotation of angle -θ around the mirror z-axis.

We obtain (r, θ) using the Gauss-Newton method by minimizing ||g|| 2 . The Implicit Function Theorem is also applied to g and provides locally a C 1function which maps (A, B) to (r, θ) and its derivation We note that all derivative calculations for one image error (Equation 4) are essentially analytical, since the "numerical part" is greatly reduced to a single 2D Gauss-Newton method.

Appendix E. Against the Exact Angular Error

This appendix shows that the exact angular error e i j = arccos(d i j .D i j ),

is not C 1 continuous when e i j = 0. Without loss of generality, the space coordinate system is changed such that d i j = ( 0 0 1 ) , and we write

x(α) y(α) z(α)


Amelia Carolina Sparavigna

An image segmentation for the measurement of microstructures in ductile cast iron

Keywords: Image processing, Segmentation of images, Size of domains, Side distribution, Ductile cast iron

The work is applying an image segmentation to a micrograph which is showing the microscopic structures of a specimen of ductile cast iron. The segmentation allows the measurement of the size of the nodules.

As explained in [START_REF] Vv | Metallography, in Wikipedia[END_REF], techniques exist to have quantitative analyses of the metallographic specimens, suitable for researches concerning metals and alloys, but also for non-metallic or composite materials. In these analyses, it is used a specimen in the form of a two-dimensional plane from a three-dimensional sample. Measurements may involve the determination of the apparent diameter of discrete particles that we can see in the specimen. To obtain such a measure, a manual approach can be used, made "with the aid of templates overlaying the microstructure" [START_REF] Vv | Metallography, in Wikipedia[END_REF], or with an image processing of the related micrographs. As stressed in Ref. 1, an adequate sampling is necessary to obtain a proper statistical basis for the measurement [START_REF] Vv | Metallography, in Wikipedia[END_REF][START_REF] Geels | Metallographic and materialographic specimen preparation, light microscopy, image analysis, and hardness testing[END_REF][START_REF] Voort | Metallography: principles and practice[END_REF][START_REF] Petzow | Metallographic etching: metallographic and ceramographic methods for reveal ing microstructure[END_REF][START_REF] Bjerregaard | Metalog guide[END_REF]. "Efforts to eliminate bias are required". Among the basic quantitative investigations, besides the measurements of the size and size distribution of particles, of their shape and of the spacing between them, we find also the measurement of the grain size in polycrystalline materials and the determination of the volume fraction of a phase or constituent. The measurement of the size and size distribution of particles is the proper task for an image segmentation method that we have recently developed and applied [START_REF] Sparavigna | Image Segmentation Applied to the Study of Micrographs of Cellular Solids[END_REF][START_REF] Sparavigna | A method for the segmentation of images based on thresholding and applied to vesicular textures[END_REF][START_REF] Sparavigna | Measuring the particles in fly ash by means of an image segmentation[END_REF][START_REF] Sparavigna | Analysis of a natural honeycomb by means of an image segmentation[END_REF][START_REF] Sparavigna | Image segmentation applied to micrographs of microcellular plastics[END_REF][START_REF] Sparavigna | Measuring the size of tubules in phloem and xylem of plants[END_REF][START_REF] Sparavigna | Image Segmentation Applied to the Analysis of Fabric Textures[END_REF]. In image processing, a segmentation is a method of partitioning an image into multiple sets of pixels, defined as super-pixels, in order to have a representation, which can be simpler than the original one or more useful to the following desired analyses [START_REF] Vv | Image segmentation[END_REF]. Specifically, the segmentation is a process of assigning a label to every pixel in an image, such that the pixels having the same label share certain characteristics [START_REF] Vv | Image segmentation[END_REF]. As a consequence, the result of a segmentation is a set of "segments", or "super-pixels", or "domains", that are covering the whole image. The segmentation can also give a set of contours, that is of "edges", extracted from the image. In this case, the segmentation gives an "edge detection". Several methods exist for segmentation, as we can appreciate from [START_REF] Vv | Image segmentation[END_REF]. Here we use a segmentation based on a thresholding of gray-scale images [START_REF] Sparavigna | Image Segmentation Applied to the Study of Micrographs of Cellular Solids[END_REF]. By selecting a suitable threshold, the image is converted into a binary (black and white) image. In several cases, this is enough for evidencing the domains among the black or the white pixels. Details of the method of segmentation are given in [START_REF] Sparavigna | Image Segmentation Applied to the Study of Micrographs of Cellular Solids[END_REF]. Let us consider the image proposed in [START_REF] Vv | Metallography, in Wikipedia[END_REF] and here reproduced in the Figure 1. It is showing a micrograph of ductile iron. This material is also known as ductile cast iron, nodular cast iron, spheroidal graphite iron, spheroidal graphite cast iron, and SG iron [START_REF] Vv | Ductile iron, in Wikipedia[END_REF][START_REF] Smith | Foundations of Materials Science and Engineering[END_REF]. It is a type of graphite-rich cast iron discovered in 1943 by Keith Millis [START_REF] Vv | Ductile iron, in Wikipedia[END_REF]. As told in [START_REF] Vv | Ductile iron, in Wikipedia[END_REF], in the Figure 1 we can see the carbon islanding effect around nodules. Actually, from the Figure 1, we can try to determine the nodules distribution, but also the areas surrounding them. Let us start from the nodules and their size distribution. As previously told, the method proposed and applied in [START_REF] Sparavigna | Image Segmentation Applied to the Study of Micrographs of Cellular Solids[END_REF][START_REF] Sparavigna | A method for the segmentation of images based on thresholding and applied to vesicular textures[END_REF][START_REF] Sparavigna | Measuring the particles in fly ash by means of an image segmentation[END_REF][START_REF] Sparavigna | Analysis of a natural honeycomb by means of an image segmentation[END_REF][START_REF] Sparavigna | Image segmentation applied to micrographs of microcellular plastics[END_REF][START_REF] Sparavigna | Measuring the size of tubules in phloem and xylem of plants[END_REF][START_REF] Sparavigna | Image Segmentation Applied to the Analysis of Fabric Textures[END_REF] is based on the thresholding of the gray tones of the image under study, to render it into a binary one. Let us remember that the values of the gray tones are integers ranging from 0 (black) to 255 (white). The histogram of the image is given in the right panel of the Figure 2. We have obtained it by means of the image processing software GIMP.

Figure 2: Screenshot of the image and its histogram obtained using GIMP.

At a first glance, the micrograph in the Figure 1 seems being an image having a trimodal histogram, that is, a histogram with three peaks. Actually, as we can see from Fig. 2, the histogram has two two evident peaks corresponding to the dark and bright areas, and between then a less pronounced peak for the gray pixels. Let us consider the peak at dark values; we can try to use two thresholds to segment the image in order to have only the pixels pertaining to this peak (for a multilevel thresholding based on the image entropy, see [START_REF] Sparavigna | Tsallis entropy in bi-level and multi-level image thresholding[END_REF]). Using GIMP we can select manually a part of the histogram, for instance, that corresponding to the blue band in the Figure 3. The selection of the gray tones in this blue stripe of the histogram are producing the binary image that we can see in the Figure 3. Let us invert the tones of the image, obtaining the Figure 4 on the left. We can compare it to the original image. We see that the nodules have their sizes slightly reduced. For this reason, we have to apply some processing before determining the size distribution. We can use the GIMP Dilate filter, which widens and enhances dark areas of the active layer or selection. We obtain the image in the Figure 5 (left panel). Using Paint, we have removed some spurious domains, having the image in Figure 5 (right panel). On this last image, we apply the segmentation [START_REF] Sparavigna | Image Segmentation Applied to the Study of Micrographs of Cellular Solids[END_REF]. The result is given in the Figure 6, which is showing the black domains as "super-pixels", each colored by a different color. After having determined the colored domains, it is easy to obtain the cross-section of each of them. The cross-section is given as the number of pixels contained in the domain. If we have the possibility to know the measure, in microns, of the image size, it is very simple to estimate the area (in microns) of the domains. Moreover, we can give the size distribution of the domains, that is, the distribution of the crosssections of nodules displayed by the micrograph of Figure 1. We obtain the distribution by counting the cross-sections according to their area (in pixels) within intervals spaced of a certain amount of pixels. Here the areas are spaced of 50 pixels, to obtain the distribution given in the Figure 7. Here we have proposed an approach based on the segmentation described in [START_REF] Sparavigna | Image Segmentation Applied to the Study of Micrographs of Cellular Solids[END_REF]. Let us conclude remarking that several software and methods exist for the image processing of micrographs for metallography. Some of them are reported in Refs. [START_REF] Maragos | Tutorial on advances in morphological image processing and analysis[END_REF][START_REF] Tewari | Application of three-dimensional digital image processing for reconstruction of microstructural volume from serial sections[END_REF][START_REF] Sternberg | Grayscale morphology[END_REF][START_REF] Stampfl | Determination of the fracture toughness by automatic image processing[END_REF].

Figure 1 :

 1 Figure 1: Light optical micrograph of polished and etched ductile iron. Etched in 3% Nital. 100X. Courtesy Michelshock, McGill University, for Wikipedia.

Figure 3 :

 3 Figure 3: Segmentation of the histogram and of the image for evidencing the nodules.

Figure 4 .

 4 Figure 4. After the segmentation and the inversion of colors, we have the image on the left, that we can compare to the image in the Figure 1.

Figure 5 :

 5 Figure 5: The image on the left is obtained by using the GIMP Dilate filter on the binary image in the Figure 4. This GIMP filter widens and enhances dark areas of the active layer or selection. Using Paint, some spurious black domains are removed (right panel).

Figure 6 :

 6 Figure 6: Result of the segmentation of the image in the Figure 1. Each nodule (super-pixel) has a different color.

Figure 7 :

 7 Figure 7: Distribution of the cross-sections of nodules. The areas are spaced in intervals of 50 pixels.

Figure 8 :

 8 Figure 8: Segmentation of the histogram and of the image for evidencing the islanding effect.