MULTI-VIEW 3D CIRCULAR TARGET RECONSTRUCTION WITH UNCERTAINTY ANALYSIS

The paper presents an algorithm for reconstruction of 3D circle from its apparition in n images. It supposes that camera poses are known up to an uncertainty. They will be considered as observations and will be refined during the reconstruction process. First, circle apparitions will be estimated in every individual image from a set of 2D points using a constrained optimization. Uncertainty of 2D points are propagated in 2D ellipse estimation and leads to covariance matrix of ellipse parameters. In 3D reconstruction process ellipse and camera pose parameters are considered as observations with known covariances. A minimal parametrization of 3D circle enables to model the projection of circle in image without any constraint. The reconstruction is performed by minimizing the length of observation residuals vector in a non linear Gauss-Helmert model. The output consists in parameters of the corresponding circle in 3D and their covariances. The results are presented on simulated data.


INTRODUCTION
Ellipses have been adopted in many computer vision applications.A 3D elliptical shape is projected as an ellipse under projective transformation.This property made detection of elliptical features easier and incited many researches in image-based object detection (Kanatani and Ohta, 2004;Martelli et al., 2010;Shuhua and Yong, 2010).Under projective projection a circle is deformed to an ellipse whenever the pose is not fronto-parallel.This property have been used for image rectification purposes (Ip and Chen, 2005;Hutter and Brewer, 2009).When 3D parameters of elliptical features and their projection in image space are known, the projection parameters can be estimated.Many researchers adopted this strategy for camera calibration (Tarel and Gagalowicz, 1995;Heikkila, 2000;Mateos, 2000;Hu and Ji, 2001).
Many man made objects exhibit circular shapes.Automatic detection of most of these features such as fiducial markers (Bergamasco et al., 2011), traffic signs (Fu and Huang, 2010;Arlicot et al., 2009), manholes (S.Ji and Shi, 2012) and roundabouts (Ravanbakhsh and Fraser, 2009) in optical images have been investigated.In contrast to 2D estimation, 3D reconstruction of circle was rarely investigated.The goal of this paper is to provide a generic approach for reconstructing circular features from their elliptical apparitions in multi-view images.The observations of our system would consist in parameters of ellipses in image space and image poses.Propagation of observations' uncertainty through the reconstruction process would provide 3D circle's uncertainty.Quan (1996) proposed a closed-form approach for conic matching and reconstruction from two views.The main drawbacks of the method consist in complexity of extension to multi-view and integration of uncertainty in the mathematical model.Moreover the provided solution is generally a 3D ellipse.Adding the circularity constraint to the solution doesn't seems to be straightforward.Mai et al. (2010) proposed a multi-view approach to estimate an ellipse in 3D space from an uncalibrated set of images.The method is based on reconstructing more than five 3D points (selected on a reference 2D ellipse) by minimizing the distances from their projections to the measured 2D ellipses on different images.The main advantages of the proposed method are in the joint analysis of multiple images (instead of stereo) and in the simultaneous estimation of pose parameters.Comparing to the previous approach, the drawback is in using representative points for 3D estimation instead of using the totality of ellipse form.Like the previous approach the solution is generally an ellipse.Lastly, the solution is not symmetric, as it depends on the choice of a reference image.Bergamasco et al. (2012) proposed a 3D parametrization of ellipses and its back-projection in image space.Then the 3D ellipse evolve in 3D space and its back-projection in multiple images are simultaneously updated.A level set function allows to fit the back-projections of the 3D ellipse to observations and define an energy.The energy is minimized by a gradient descent method.Comparing to the previous multi-view method, this approach ensure the consistency of the back-projection in all ellipses' points in any image.
The former method is the nearest method to ours except in formulation for circle and taking into account the uncertainty.Uncertainty propagation in terms of projective geometries has been developed by Förstner (2005).We propose to adapt these concepts of error propagation in our context of 3D circle reconstruction under uncertain views.Uncertainty analysis and propagation may be carried out using the Gauss-Helmert model (McGlone et al., 2004), which enables simultaneous estimation of unknowns and adjustments of observations while providing uncertainty propagation in terms of covariance matrices.

Proposed approach
To tackle the simultaneous 3D circle reconstruction problem and multi-view 3D pose correction problem with error propagation, we propose a stratified approach with 2 steps : 2. These estimated 2D ellipses with uncertainties are then taken as observations to be adjusted in the subsequent 3D reconstruction step which estimates a 3D circle and adjusts the poses of all images, carrying out uncertainty propagation on all these elements (section 3).
Section 4 and 5 provide some results, discussion and conclusion on the proposed method.

2D geometry concepts and notations
A conic is a 2D shape that may be described by the following implicit equation : This may be rewritten in homogeneous coordinates with a 3 × 3 symmetric matrix E (Hartley and Zisserman, 2004) : Dual to this point based representation, a conic may also be represented as a curve that is tangent to a set of lines.Using homogeneous line coordinates ( a, b, c ), the tangency relationship may be written as , where E * is the line conic dual to the point conic E. If E is full-rank, then we have the following relationship, as E and E * are symmetric homogeneous quantities : This dual representation will not be used to characterize tangent lines but for the simplicity expressing a projected 3D dual quadric as a 2D dual conic (cf.section 3.2).
Finally, ellipses are the subset of conics for which ac−b 2 > 0. As the symmetric matrix E is a homogeneous quantity defined up to scale, encoding an ellipse may thus be scaled to fulfill ac − b 2 = 1.

Ellipse fitting and error propagation
Estimation of 2D ellipse consists in resolution of the following constrained equation system : where: x : (a, b, c, d, e, f ) T : unknowns vector l : (x1, y1, ..., xn, yn) T : observations vector We applied the Gauss-Helmert model to resolve the system by minimizing length of observations residual vector.The model enables to take into account covariance matrix of observations to weight the equations and provides the covariance matrix of the parameters (Vaníček and Krakiwsky, 1986).As the constraint function fc is not linear, both equations were linearized to yield : x 0 , l 0 : Initial values δ, r : residuals of unknowns and observations The variation function for finding the least-squares solution is : The unknowns correction vector δ and observations correction vector r are obtained by minimizing Φ (Förstner, 2005). where: Cr : Covariance matrix of observations The Gauss-Newton iterative method were employed to solve the problem by applying corrections to the observations r and unknowns δ iteratively.We applied the method presented by Fitzgibbon et al. (1999) in order to obtain an initial solution x 0 .
Finlly the covariance matrix of parameters C δ is obtained from inverted normal equation matrix:

Point Set Normalization
To prevent numerical issues, the point set (Xi) = (xi, yi) is converted through a linear transform to a normalized point set (X i ) = (x i , y i ), prior to the Ellipse fitting and error propagation.The inverse of this linear transform, denoted by the homogeneous denormalization matrix M den verifies Xi = M den X i and may be written as : where (sx, sy) is the normalization scale and (mx, my) the normalization center.
The results of ellipse fitting is the ellipse E and covariance Σ E estimated on the normalized point set.

Ellipse Dualization
The estimated dual conic E * is proportional to the inverse of E , and may thus be chosen as the comatrix of E (equation 3).Since ac − b 2 = 1, we get : Uncertainty propagation may thus be carried out with the following approximation : with

Dual Ellipse Denormalization
The dual conic E * may finally be derived from the dual conic E * of the normalized point set using the following identity, expressing a dual conic modification under the point transform M den : This amounts to a linear transform on the vector representation vec5(E * ) of E * , which may be used to compute an exact uncertainty propagation on its covariance matrix. with Chaining all these steps together, we have thus proposed a method that estimates the uncertain dual ellipse quantities (vec5(E * ), Σ vec 5 (E * ) ) from a set of ellipse contour points in the 2D image.

Dual Quadric of a 3D Circle
Projective geometry may also be used to model (point) quadrics and dual (plane) quadrics (Hartley and Zisserman, 2004).A 3D circle may not be modeled using a quadric Q, defined by homogeneous points X such that X T QX = 0.It may however be defined using a dual quadric Q * , defined by homogeneous planes Π such that Π T Q * Π = 0.These planes tangent to a 3D circle are the planes tangent to its rim, including its supporting plane.By denoting C the center point of the 3D circle and N a unit vector normal to its supporting plane scaled by the circle radius ρ (i.e.ρ 2 = N 2 ), the dual quadric Q * is singular (rank 3) and may be written as : where [N ]× is the matrix encoding the vector product [N ]×X = N × X.This proposed 6D parameterization (C, N ) is a minimal parametrization of 3D circles that is both unconstrained and unambiguous except for the sign of N .
Proof By squashing the quadric Qt = diag(1, 1, t, −ρ 2 ) from the origin-centered sphere of radius ρ at t = 1 to a disc of radius ρ in the z = 0 plane when t → ∞, we get its dual quadric Q * ∞ : ] is a rotation matrix with N = ρW .Given the transformation rule of dual quadrics : which yields equation 23.

Dual Conic of a Projected 3D Circle
A dual quadric Q * is imaged as a dual conic E * ∼ = P Q * P T by a projection P .The image projection P is usually decomposed as KR I3 −S with S denoting the projection center, R the projection rotation and K the intrinsic matrix .Using equation 23, this yields, up to a scale factor and defining M = KR, the dual conic of a projected 3D circle (C, N ) : By splitting row-wise M = [ M 1 M 2 M 3 ] T , equation 26 may be rewritten in terms of the elements of E * :

Resolution of the equation system
Each image observation of a 2D ellipse introduces a set of equations translating that the 2D ellipse projected from the 3D circle E * (equation 26) should agree with the 2D ellipse estimated from ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-3, 2014ISPRS Technical Commission III Symposium, 5 -7 September 2014, Zurich, Switzerland This contribution has been peer-reviewed.The double-blind peer-review was conducted on the basis of the full paper.doi:10.5194/isprsannals-II-3-143-2014 image contours E * obs (section 2).This induces a set of 5 equations and no constraints, thanks to the proposed unconstrained parametrization.Namely, E * obs must verify, up to a scalar factor, equation 26, which yields 6 equations (E * obs being a 3 × 3 symmetric matrix) minus one for the unknown homogeneous scale factor λ. Given that E * obs,33 = 1, this writes : E * obs is considered here as an observation with the uncertainty Σ E * obs .E * may be derived from the observations M, S and the 3D circle unknowns C, N .
Each 2D ellipse observation thus provides a 5-dimensional observation vector.The pose parameters of the corresponding image add 6 other observations (3 coordinates of the projection center and 3 rotations).We suppose that covariances matrix of pose parameters are known.Intrinsic parameters could be considered as observation in the same way but in our implementation we considered them as fix values.The observation vector provided by each image apparition become: The 3D circle to be estimated is encoded in the 6-vector X = (C, N ).In this context, the 5-dimensional equation 28 may be rewritten in terms of the observations L and the unknowns X : F and E * feature the following derivatives, introducing the 3 × 3 singular matrices A ij = (MiM T j + MjM T i ) and the canonical vectors δ i which are 0 everywhere but 1 at the i th element : where From these equations, it is straightforward to derive the Jacobian matrices ∂F ∂X = [ ∂F ∂C ∂F ∂N ] and ∂F ∂L .The system is resolved using unconstrained Gauss-Helmert model who can be obtained by replacing D = 0 in equation 9.The nonlinear system is again resolved by Gauss-Newton iterative method.The output of this step consists in adjusted values of unknowns ( Ĉ, N ) and their uncertainties (C Ĉ , C N ) together with adjusted observations l.

Data simulation
We set up two scenarios of three perspective images of size 1000× 2000 pixels and a focal length of 1400 pixels.A circular disk of 80cm of diameter was placed at a mean distance of 10 meters from the images (cf.Fig. 1).The first scenario was composed of images Set1 : (I − III − IV ) and the second one images Set2 : (II − III − IV ).The approximate distances between cameras was 4m.The difference between the two sets is that in the former the camera centers are almost on the same plane whereas in the latter the camera II was at a distance of 3m from the horizontal plane containing the two other cameras.The 3D circular disk was sampled regularly at ten 3D points and projected in each camera.The observations consisted in : • 2D noisy coordinates of projected sampled points in the images.
• 3D noisy projection centers of the cameras.
• Rotation matrices of the cameras.

Simulation results
For each scenario we added a Gaussian noise to the projection centers of cameras.Then we run both steps of ellipse estimation and 3D reconstruction 20 times for each noise level.The output of each iteration consists in the estimated values of the parameters and their covariance matrix.Using covariance matrices, 3D error ellipsoids (99%) were computed for the center and normal vector of the circle.Fig. 2 depicts on of these iterations.Knowing the reference parameters, the difference between the obtained parameters and the ground trough can be computed.
Ideally the probability that the computed parameters match the reference parameters could be computed using the estimated covariances.However in this paper in order to show at the same time how the estimated parameters vary as a function of the applied noise we show in the same diagrams the norm of difference vector and the length of largest error ellipsoid axis for both center and normal vector for both scenarios.Fig. 3 shows the results for Both diagrams demonstrate that generally the errors are larger in the first scenario than the second for both center and normal vector.This is natural since non aligned cameras provide better 3D intersections.This fact is also confirmed by the estimated error ellipsoids since that of second scenario is smaller than the second one.This is also visible on Fig 2 for one iteration where a Gaussian error of σ = 1cm were added to the projection center of cameras.The effective difference between the estimated parameter and the reference one is always smaller than the largest axis of error ellipsoid (99%).This means that the estimated precision of the network is confirmed by the real precision measured on control points.
On can understand from this experiment that the second image network is much less sensible to pose errors.It can theoretically lead to 10cm of error in 3D for a pose error of 4cm whereas the first image network for the same amount of pose error can lead to 60cm of error in 3D.By a simple repositioning of cameras on can improve the precision of the 3D measurement system.On can also understand that the improvement (of using scenario 2 vs scenario 1) is much more significant for estimation of normal vector than that center of the circle.

CONCLUSIONS AND PERSPECTIVES
We proposed a new approach for reconstructing circular targets from uncertain multiple-views.Our main contributions are in estimating the uncertainty of 2D ellipse and in formalizing perspective projection of circle with minimum number of parameters  It can be applied for reconstruction of any circular feature using an unlimited number of images with uncertain calibration.It enables also to estimate the precision and stability of a photogrammetric network for measuring circular features.Therefore position of the cameras can be optimized (through an analysissynthesis schema) in order to improve the 3D precision without multiplying the number of the cameras.
The method should be evaluated on more simulated scenarios.Namely by studying its behavior facing more noises on initial 2D points (used for ellipse fitting) and camera rotation matrix but also with more images in complicated camera configurations.Then we aim at evaluating its performance on real data in reconstruction of circular traffic signs (in the same framework developed for reconstruction of polygonal road signs (Soheilian et al., 2013)) and also in detection and reconstruction of circular photogrammetric calibration targets.
More theoretically we aim at adapting the method to handle multiple circular features reconstruction at the same time.In more long term we will integrate the developed method in a bundle adjustment chain enabling use of circular features together with their uncertainties as tie and control points.

Figure 1 :
Figure 1: Configuration of data simulation

Figure 2 :
Figure 2: Ellipsoids of errors (99%) corresponding to the center and normal for both scenarios.Noise of projection center = σpc = 1cm normal vector and Fig. 4 depicts that for center.
error applied to the projection center (m) Scenario 1: difference Scenario 1: Length of largest error ellipsoid axis Scenario 2: difference Scenario 2: Length of largest error ellipsoid axis

Figure 3 :
Figure 3: Error induced by random errors of projection center to normal vector.