Direct Least Squares Fitting of Ellipses

Andrew W. Fitzgibbon Maurizio Pilu Robert B. Fisher Department of Articial Intelligence The University of Edinburgh 5 Forrest Hill, Edinburgh EH1 2QL SCOTLAND andrewfg@[Link]

Abstract
This work presents a new efcient method for tting ellipses to scattered data. Previous algorithms either tted general conics or were computationally expensive. By minimizing the algebraic distance subject to the constraint 4 2 1 the new method incorporates the ellipticity constraint into the normalization factor. The new method combines several advantages: (i) It is ellipse-specic so that even bad data will always return an ellipse; (ii) It can be solved naturally by a generalized eigensystem and (iii) it is extremely robust, efcient and easy to implement. We compare the proposed method to other approaches and show its robustness on several examples in which other non-ellipse-specic approaches would fail or require computationally expensive iterative renements. Source code for the algorithm is supplied and a demonstration is available on !#"$&%('!)01%324'%3)657%98@A4BA)#8CDFE

GIHH

QPDFRA24BA"'!24BA"7%(SB

Introduction

conic and rejecting non-elliptical ts. These latter methods are cheap and perform well if the data belong to a precisely elliptical arc with little occlusion but suffer from the major shortcoming that under less ideal conditions non-strictly elliptical data, moderate occlusion or noise they often yield unbounded ts to hyperbolae. In a situation where ellipses are specically desired, such ts must be rejected as useless. A number of iterative renement procedures [18, 8, 13] alleviate this problem, but do not eliminate it. In addition, these techniques often increase the computational burden unacceptably. This paper introduces a new tting method that combines the advantages of: Ellipse-specicity yielding useful results under all noise and occlusion conditions; Afne invariance; and Robustness to noise and occlusion. After a description of previous algebraic tting methods, in Section 3 we describe the method and provide a theoretical analysis of the uniqueness of the elliptical solution. Section 4 contains experimental results, notably to highlight noise resilience, invariance properties and behaviour for non-elliptical data. We conclude by presenting some possible extensions.

The tting of primitive models to image data is a basic task in pattern recognition and computer vision, allowing reduction and simplication of the data to the benet of higher level processing stages. One of the most commonly used models is the ellipse which, being the perspective projection of the circle, is of great importance for many industrial applications. Despite its importance, however, there has been until now no computationally efcient ellipse-specic tting algorithm [15, 6]. In this paper we introduce a new method of tting ellipses, rather than general conics, to segmented data. As we shall see in the next section, current methods are either computationally expensive Hough transform-based approaches, or perform ellipse tting by least-squares tting to a general

Previous Methods and their Limitations

The literature on ellipse tting divides into two general techniques: clustering and least-squares tting. Clustering methods are based on mapping sets of points to the parameter space, such as the Hough transform [10, 22] and accumulation methods [17]. These Hough-like techniques have some great advantages, notably high robustness to occlusion and no requirement for pre-segmentation, but they suffer from the great shortcomings of high computational complexity and non-uniqueness of solutions, which can render them unsuitable for real applications. Particularly when curves have been pre-segmented, their computational cost is signicant.

Least-squares techniques center on nding the set of parameters that minimize some distance measure between the data points and the ellipse. In this section we briey present the most cited works in ellipse tting and its closely related problem, conic tting. It will be shown that the direct specic least-square tting of ellipses has, up to now, not been solved.

A B

2.1

Problem statement

   2     2 !  #"  #$ 0 (1) 2%  2 3' 2 4 (10 where &% '"$)(10 and 5 ; 768 is called the algebraic distance of a point  1  . ;  0. The tting of a general conic to the conic
may be approached [7] by minimizing the sum of squared algebraic distances

Before reviewing the literature on general conic tting, we will introduce a statement of the problem that allows us to unify several approaches under the umbrella of constrained least squares. Let us represent a general conic by an implicit second order polynomial:

Figure 1. Hand-drawn data. The linetype / algorithm correspondences are Bookstein: dotted; Taubin: dash-dot; Ellipse-specic: solid.

@ 9 768 2 (2) 65A 1 6 of the curve to the B data points . In order to avoid the trivial solution DC 6 , and recognizing that any multiple represents the same conic, the parameter of a solution vector is constrained in some way. Many of the published
algorithms differ only in the form of constraint applied to the parameters:

E E E E E E

2 Many authors [11] suggest 1. Bolles and Fischler [2] suggest 1. 1. Rosin [15] and Gander [6] impose Rosin [15] also investigates 1. 2 1. 1 2 2 Bookstein [3] proposes 2 Agin [1] and Taubin [20] use the data-dependent quad2 ratic constraint 1 where is the Jacobian ; 1 ; .

F F $

F G F B Q G P R P P I% H H 5 9 1( 0

Note that these constraints are all either linear, of the form 1 or quadratic, constraining 1 where is a 6 6 constraint matrix. A number of papers have concerned themselves with the specic problem of recovering ellipses rather than general conics. The above methods do not restrict the tting to be an ellipse, in the sense that they can return an hyperbola or a parabola even given elliptical data. Porrill [13] and Ellis et al. [4] use Booksteins method to initialize a Kalman lter. The Kalman lter iteratively minimizes the gradient distance in order to gather new image

S 

0UT

evidence and to reject non-ellipse ts by testing the discrim0 at each iteration. Porrill also gives nice inant 2 4 examples of the condence envelopes of the ttings. Although these methods transform the disadvantage of having a non-specic ellipse tting method into an asset by using the ellipse constraint to check whether new data has to be included or to assess the quality of the t, the methods require many iterations in the presence of very bad data, and may fail to converge in extreme cases. Rosin [15] initializes a Kalman Filter with a circle, and tightens the corresponding initial covariance until the returned conic is an ellipse. In [16] he restates that ellipsespecic tting is a non-linear problem and that iterative methods must be employed. He also [15] analyses the pro and cons of two commonly used normalizations, 1 and 1 and shows that the former biases the tting to have smaller eccentricity, therefore increasing the probability of returning an ellipse, at the cost of losing transformational invariance. Haralick [7, 11.10.7] takes a different approach. Effectively, he guarantees that the conic is an ellipse by replacing the coefcients with new expressions 2 2 2 2 so that the discriminant 2 4 be 2 2 comes 4 which is guaranteed negative. Minimization over the space then yields an ellipse. His algorithm is again iterative, and an initial estimate is provided by a method of moments. Keren et al. [9] apply a similar technique to Haralicks in the context of the tting of bounded quartic curves. Again, their algorithm is iterative.

XW

ba `  f a `c  d ce e g c f f "$ a `dc e

Direct ellipse-specic tting

In order to t ellipses specically while retaining the efciency of solution of the linear least-squares problem (2), we would like to constrain the parameter vector so that the conic that it represents is forced to be an ellipse. The appropriate constraint is well known, namely that the discriminant 2 4 be negative. However, this constrained problem is difcult to solve in general as the Kuhn-Tucker conditions [14] do not guarantee a solution. In fact, we have not been able to locate any reference regarding the minimization of a quadratic form subject to such a nonconvex inequality. Although imposition of this inequality constraint is difcult in general, in this case we have the freedom to arbitrarily scale the parameters so we may simply incorporate the scaling into the constraint and impose the equality constraint 4 2 1. This is a quadratic constraint which may be expressed in 1 as the matrix form

10

10

10

10

Figure 2. Average geometric distance error as a function of increasing noise level. The errorbars are at ( 1 ) . The pictures along the noise axis indicate visually the corresponding noise level. Encoding is Bookstein: dotted; Taubin: dash-dot; New: solid.

0 T
0 0 2 0 0 0

0 1 0 0 0 0

2 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

(3)

3.1

Solution of the quadratically constrained minimization

Following Bookstein [3], the constrained tting problem is: Minimize

F F 2, subject to 0 T


where the design matrix is dened as the  6 matrix   . Introducing the Lagrange multiplier 1 2 and differentiating we arrive at the system of simultaneous equations1 2

% RG ( 0

(4)

2

T 0 T 

0 1

solves (5). As in general there Finally, setting may be up to 6 real solutions, the solution is chosen that 21 yields the lowest residual  . We note that the solution of the eigensystem (6) gives 6  eigenvalue-eigenvector pairs  . Each of these pairs gives rise to a local minimum if the term under the square  root in (8) is positive. In general, is positive denite, so  "  the denominator is positive for all . Therefore the 43 square root exists if  0, so any solutions to (5) must have positive generalized eigenvalues. A Cholesky decom1 position 65 2 reveals that the signs of the eigenvalues of (6) are the same as those of 57 1 845' 1, and therefore, by Sylvesters law of Inertia [21] the same as those of 8 . The constraint matrix (3) for the ellipse-specic problem has only one positive eigenvalue, and therefore the elliptical solution is unique [12].
 0

6

76 6 06 6 6 6d 68 6 6

06

(5)

3.2
(6) (7)

Observation

This may be rewritten as the system

7 0 T


Solving (4) is equivalent to minimizing the Rayleigh quo@9 9 tient , which in this case is the cost function

where is the scatter matrix . This system is readily solved by considering the generalized eigenvectors of (6). If     solves (6) then so does  for any and from     1 giving (7) we can nd the value of as 2

6 6

 

6 6 6 6 06 T 6 6 1 06 T 6 06 6
 !   "

@ 9 5 6 2 2 65A 1

(8)

1 Note that the method of Lagrange multipliers is not valid when the $# gradient of& the constraint function becomes zero. In (4) this means % '# 0, but then 0 so the constraint is violated and there is no solution.

Now, as the discriminant is inversely proportional to the product of the radii, this is analogous to multiplying the cost function by the area of the ellipse. Interestingly, this is the technique adopted by Solina and Bajcsy for tting superquadrics to range data [19], in order to bias the t towards smaller volume. This provides an intuitive feeling for the low-curvature bias of the method.

Experimental Results
% x,y are lists of coordinates function a = fit ellipse(x,y) % Build design matrix D = x.*x x.*y y.*y x y ones(size(x)) ; % Build scatter matrix S = D*D; % Build 6x6 constraint matrix C(6,6)=0; C(1,3)=2; C(2,2)=-1; C(3,1)=2; % Solve eigensystem gevec, geval = eig(inv(S)*C); % Find the positive eigenvalue PosR, PosC = find(geval 0 & isinf(geval)); % Extract eigenvector corresponding to positive eigenvalue a = gevec(:,PosC);

In this section we present experimental results that compare the ellipse-specic solution to previous methods in terms of quality and robustness. We include both quantitative and qualitative results in order to allow other researchers to evaluate the utility of the ellipse-specic algorithm with respect to the others cited. Fitzgibbon [5] provides further theoretical and quantitative results for a wide range of conic-tting algorithms.

4.1

Ellipse-specicity
Figure 4. Complete 6-line Matlab implementation of the proposed algorithm.

Despite the theoretical proof of the algorithms ellipsespecicity, it is instructive to observe its performance on some real data, of which Figure 1 provides some examples with hand-drawn datasets. The results of our method are superimposed on those of Bookstein and Gander. Dataset A is almost elliptical and indistinguishable ts were produced. The other sets exhibit varying degrees of non-ellipticity, and illustrate the potential use of ellipse-specicity for coarse data bounding of nonelliptical data.

4.3

Afne invariance

4.2

Noise sensitivity

The rst noise experiment measures the average geometric distance error for each of the algorithms over 100 runs. In order to verify that the ellipses returned by the new algorithm are reasonable approximations to the minimum geometric distance ellipse, non-elliptical ts returned by the Bookstein and Taubin algorithms were ignored. It can be seen that our algorithm produces a closer ellipse on average than Booksteins for medium noise, but that Taubinswhen it returns an ellipseproduces the smallest geometric distance error. We note however that all results are within each others 1 ) error bars over the 100 runs, meaning that the variations within runs are greater than the difference between the algorithms across runs. The second experiment, illustrated in Figure 3, is perhaps more important (although we have not seen it in related papers) and is concerned with assessing the stability of the tting with respect to different realizations of noise with the same variance. It is very desirable that the algorithm performance be affected only by the noise level, and not by a particular realization of the noise. Figure 3 shows ten different runs in which a different noise population with same variance ( ) 0 1) was generated and results for each of the three methods is displayed. In this and similar experiments we found that the stability of the method is noteworthy. Ganders algorithm shows a greater variation in results and Taubins, while improving on Ganders, remains less stable than the proposed algorithm.

The constraint 4 2 1 not only constrains the tted conics to be ellipses but it is also covariant with afne transformations of the data points. In an experimental determination, we applied random nonsingular afne transforms to the points of a data set, applied the tting algorithm, inverse-transformed the ellipse parameters by the original transform and compared the recovered parameters to the expected ones. The difference was zero up to machine precision.

Conclusions

This paper has presented a new method for direct least square tting of ellipses. We believe this to be the rst noniterative ellipse-specic algorithm. Previous conic tting methods rely (when applied to ellipse tting) either on the presence of good data or on computationally expensive iterative updates of the parameters. We have theoretically demonstrated that our method uniquely yields elliptical solutions that, under the normalization 4 2 1, minimize the sum of squared algebraic distances from the points to the ellipse. Experimental results illustrate the advantages conferred by ellipse specicity in terms of occlusion and noise sensitivity. The stability properties widen the scope of application of the algorithm from ellipse tting to cases where the data are not strictly elliptical but need to be minimally represented by an elliptical blob. In our view, the method presented here offers the best tradeoff between speed and accuracy for ellipse ttingits simplicity is demonstrated by the inclusion in Figure 4 of a complete 6-line implementation in MATLAB. In cases where more accurate results are required, this algorithm provides an excellent initial estimate. We note also that the algorithm

Sigma=0.10

Sigma=0.10

Sigma=0.10

Sigma=0.10

Sigma=0.10

Sigma=0.10

Sigma=0.10

Sigma=0.10

Sigma=0.10

Sigma=0.10

Figure 3. Stability experiments for different runs with same noise variance. Top row: proposed method; Mid Row: Ganders Method; Bottom Row: Taubins method

can be trivially converted to a hyperbola-specic tter, and a variation may be used to t parabolae. The algorithm is however biased towards ellipses of low eccentricity, and future work includes the incorporation of the algorithm into a bias-correction algorithm based on that of Kanatani [8]. In a similar vein, a theoretical analysis of the noise performance of the methods using eigensystem perturbation theory is under investigation.

6 Acknowledgements
The second author is partially sponsored by SGSTHOMSON Microelectronics UK. This work was partially funded by UK EPSRC Grant GR/H/86905.

References
[1] G. J. Agin. Representation and description of curved objects. PhD thesis, AIM-173, Stanford AI Lab, 1972. [2] R. C. Bolles and M. A. Fischler. A RANSAC-based approach to model tting and its application to nding cylinders in range data. In Proceedings, International Joint Conference on AI, pages 637643, 1981. [3] F. Bookstein. Fitting conic sections to scattered data. Computer Graphics and Image Processing, (9):5671, 1979. [4] T. Ellis, A. Abbood, and B. Brillault. Ellipse detection and matching with uncertainty. Image and Vision Computing, 10(2):271276, 1992. [5] A. Fitzgibbon and R. Fisher. A buyers guide to conic tting. In Proceedings of British Machine Vision Conference, Birmingam, 1995. [6] W. Gander, G. Golub, and R. Strebel. Least-square tting of circles and ellipses. BIT, (43):558578, 1994. [7] R. M. Haralick and L. G. Shapiro. Computer and Robot Vision, volume 1. Addison-Wesley, 1993. [8] K. Kanatani. Statistical bias of conic tting and renormalization. IEEE T-PAMI, 16(3):320326, 1994.

[9] D. Keren, D. Cooper, and J. Subrahmonia. Describing complicated objects by implicit polynomials. IEEE T-PAMI, 16(1):3853, 1994. [10] V. Leavers. Shape Detection in Computer Vision Using the Hough Transform. Springer-Verlag, 1992. [11] K. A. Paton. Conic sections in automatic chromosome analysis. In B. Meltzer and D. Michie, editors, Machine Intelligence 5, pages 411434. Edinburgh University Press, Edinburgh, 1969. [12] M. Pilu. Part-based Grouping and Recogntion: A ModelGuided Approach. PhD Thesis, Department of Articial Intelligence, University of Edinburgh, Scotland, 1996. Jun, 1996. [13] J. Porrill. Fitting ellipses and predicting comdence envelopes using a bias corrected Kalman lter. Image and Vision Computing, 8(1), Feb. 1990. [14] S. S. Rao. Optimization: Theory and Applications. Wiley Eastern, 2 edition, 1984. [15] P. Rosin. A note on the least square tting of ellipses. Pattern Recognition Letters, (14):799808, Oct. 1993. [16] P. Rosin and G. West. Segmenting curves into lines and arcs. In Proceedings of the Third International Conference on Computer Vision, pages 7478, Osaka, Japan, Dec. 1990. [17] P. L. Rosin. Ellipse tting by accumulating ve-point ts. Pattern Recognition Letters, (14):661699, Aug. 1993. [18] P. Sampson. Fitting conic sections to very scattered data: An iterative renement of the Bookstein algorithm. Computer Graphics and Image Processing, (18):97108, 1982. [19] F. Solina and R. Bajcsy. Recovery of parametric models from range images: The case of superquadrics with global deformations. IEEE PAMI, 12(2), Feb. 1990. [20] G. Taubin. Estimation of planar curves, surfaces and nonplanar space curves dened by implicit equations , with applications to edge and range image segmentation. IEEE PAMI, 13(11):11151138, Nov. 1991. [21] J. H. Wilkinson. The algebraic eigenvalueproblem. Claredon Press, Oxford, England, 1965. [22] H. K. Yuen, J. Illingworth, and J. Kittler. Detecting partially occluded ellipses using the Hough transform. Image and Vision Computing, 7(1):3137, 1989.