Direct Least Square Fitting of Ellipses

Andrew FitzgibbonM. PiluRobert B. Fisher

article1999TPAMI2,764 citations

Proposes the first non-iterative least-squares method for fitting ellipses to scattered 2D data by solving a quadratically constrained generalized eigensystem that mathematically guarantees an elliptical solution even under severe noise and occlusion.

Listen

Fitting ellipses to scattered data points supports key tasks in computer vision and industrial inspection, where ellipses represent perspective projections of circles. Prior approaches either solved for general conics and often returned hyperbolas under noise or occlusion, or relied on slow iterative refinement to enforce ellipticity.

The article set out to create a direct, non-iterative least-squares method that guarantees an elliptical solution while remaining computationally cheap and invariant to Euclidean transformations.

The authors minimize the sum of squared algebraic distances subject to the quadratic constraint 4ac - b² = 1. This constraint is incorporated into a generalized eigenvalue problem whose solution yields exactly one ellipse. Experiments compared the new method against Bookstein, Gander, and Taubin algorithms on synthetic arcs with increasing Gaussian noise and on hand-drawn non-elliptical data.

The method always returns an ellipse even from poor data, exhibits graceful degradation under rising noise levels up to 17 percent of data spread, and produces more stable fits across repeated noise realizations than the baselines. It matches or exceeds the geometric accuracy of prior techniques when those techniques happen to return ellipses, runs in a few lines of code, and remains unchanged under rotation or translation of the input points. A side effect is a bias toward low-eccentricity ellipses.

These properties make the algorithm immediately usable for real-time ellipse detection and for supplying reliable initial estimates to more expensive iterative refiners. It widens applicability to cases where data are only approximately elliptical yet an elliptical summary is still required.

The authors recommend the method as a fast default or initializer and note that future bias-correction steps, already explored in related work, could remove the low-eccentricity tendency. The approach can be trivially adapted to fit hyperbolas or parabolas instead.

The main limitations are the low-eccentricity bias and reliance on algebraic rather than geometric distance; however, the theoretical uniqueness proof and consistent experimental behavior across dozens of test cases give high confidence in the core claims.

University of Edinburgh demonstrations html http Direct Least Square Fitting Ellipse Specific Generalized Eigenvalue

Fitzgibbon et al (1999).pdf

No sufficiently relevant recommendations were found.

Cover for Direct Least Square Fitting of Ellipses

Abstract

This work presents a new efficient method for fitting ellipses to scattered data. Previous algorithms either fitted general conics or were computationally expensive. By minimizing the algebraic distance subject to the constraint 4ac − b2 = 1 the new method incorporates the ellipticity constraint into the normalization factor. The new method combines several advantages: (i) It is ellipse-specific 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, efficient and easy to implement. We compare the proposed method to other approaches and show its robustness on several examples in which other non-ellipse-specific approaches would fail or require computationally expensive iterative refinements. Source code for the algorithm is supplied and a demonstration is available on http://www.dai.ed.ac.uk/groups/mvu/ellipse-demo.html

Topics

demonstrations html http source codes adaptive algorithms optimization problems

Keywords

Direct Least Square Fitting Ellipse Specific Generalized Eigenvalue Algebraic Distance Constrained Least Square Robustness Noise Sensitivity Euclidean Invariance Iterative Refinement Quadratic Constraint

1 Introduction
2 Previous Methods and their Limitations
2.1 Problem statement
2.2 General conic #0Ctting
2.3 Towards ellipse-speci#0Cc #0Ctting
3 Direct ellipse-speci#0Cc #0Ctting
3.1 Solution of the quadratically constrained minimization
3.2 Analysis of the constraint 4 , 2 = 1
3.3 Remark
4 Experimental Results
4.1 Ellipse-speci#0Ccity
4.2 Noise sensitivity
4.3 Parabolic #0Ct
4.4 Euclidean transformation invariance
5 Conclusions
6 Acknowledgements
References
Appendix

Knowls

Knowl 1 — Conic Section Fitting via Levenberg-Marquardt Optimization
model/method
fitConic fits two-dimensional point coordinates $(x_i, y_i)_{i=1}^N$ to an ellipse, hyperbola, or parabola, supporting arbitrary planar rotation.

Fitting proceeds through the following pathways:
- Ellipses and Hyperbolas: If initial parameter estimates $\mathbf{p}_0$ (given either as six algebraic coefficients $[A, B, C, D, E, F]$ or five geometric parameters $[h, v, a, b, \theta]$ ) are not provided, starting values are obtained via bootstrap estimators (bootEllipse for ellipses or bootHyperbola for hyperbolas). Parameters are iteratively refined via a Levenberg-Marquardt algorithm (LMA) to minimize geometric orthogonal residuals. Optimization is governed by initial damping $\lambda_0$ (LambdaIni, default $1$ ), parameter step tolerance $\epsilon_P$ (epsilonP, default $10^{-6}$ ), function tolerance $\epsilon_F$ (epsilonF, default $10^{-6}$ ), and maximum iterations $N_{\max}$ (IterMAX, default $20000$ ).
- Parabolas: When fitting a parabola, fitConic dispatches directly to fitParabola, which performs a global angular search and derotated polynomial fit without requiring initial parameter guesses.
- Outputs: The function returns the algebraic quadratic coefficients $\mathbf{p} = [A, B, C, D, E, F]$ , root-sum-square (RSS) residual error, iteration count at convergence, and an exit code ( $1$ for ellipse, $2$ for hyperbola, $3$ for parabola, or negative/non-positive flags for degenerate fits like straight lines).
Knowl 2 — Algebraic and Geometric Parameterizations of Conic Sections
definition
In Cartesian coordinates $(x, y)$ , conic sections are parameterized in two primary equivalent forms:
1. General Quadratic (Algebraic) Form: $A x^2 + B x y + C y^2 + D x + E y + F = 0$ represented by the coefficient vector $\mathbf{p}_A = [A, B, C, D, E, F]$ .
2. Center-Axis-Angle (Geometric) Form: Represented by $\mathbf{p}_G = [h, v, a, b, \theta]$ , where $(h, v) \in \mathbb{R}^2$ is the center, $a > 0$ is the semi-major axis, $b > 0$ is the semi-minor axis, and $\theta$ is the rotation angle: $\frac{((x - h)\cos\theta + (y - v)\sin\theta)^2}{a^2} + s \cdot \frac{((x - h)\sin\theta - (y - v)\cos\theta)^2}{b^2} = 1$ where $s = +1$ for an ellipse and $s = -1$ for a hyperbola. The full axis lengths are $2a$ and $2b$ .
The conic section type is categorized by the discriminant $\Delta = B^2 - 4AC$ :
- $\Delta < 0$ : Ellipse
- $\Delta = 0$ : Parabola
- $\Delta > 0$ : Hyperbola
Conversion functions (AtoG, GtoA) apply a tolerance threshold $\text{tol}$ (default $10^{-6}$ ); if $|B| < \text{tol}$ , $B$ is set to zero (no rotation), and if $|B^2 - 4AC| < \text{tol}$ , the curve is identified as a parabola.
Knowl 3 — Derotation Angle Formula for Quadratic Conics
equation
For a general quadratic conic section defined by $A x^2 + B x y + C y^2 + D x + E y + F = 0$ with $A, B, C, D, E, F \in \mathbb{R}$ and $B \neq 0$ , the rotation angle $\theta$ required to eliminate the cross-term $B x y$ (derotate the conic so that its axes align parallel to the Cartesian coordinate axes) satisfies: $\cot(2\theta) = \frac{A - C}{B} \implies \tan(2\theta) = \frac{B}{A - C}$ Applying this rotation yields a transformed parameter set $[A', 0, C', D', E', F']$ orthogonal to the Cartesian axes.
Knowl 4 — Rotated Parabola Fitting Algorithm
algorithm
The fitParabola function estimates the parameters of an arbitrarily rotated parabola from a set of 2D data points without requiring an initial parameter guess.
Input: Point coordinates $X = (x_i, y_i)_{i=1}^N$ , angular search bounds $[ heta_{\min}, \theta_{\max}]$ (default $[-\pi/2, \pi/2]$ )
Output: Vertex coordinates $(x_v, y_v)$ , rotation angle $\theta^*$ , algebraic coefficients $[A, B, C, D, E, F]$ , polynomial coefficients $[c_2, c_1, c_0]$ , cost
Find optimal angle $\theta^* = \arg\min_{\theta \in [\theta_{\min}, \theta_{\max}]} \text{cost}(\theta, X)$ using 1D optimization (RANSAC-style search)
Rotate input coordinates by $-\theta^*$ to make the parabola vertical:
     $x'_i = x_i \cos(-\theta^*) - y_i \sin(-\theta^*)$
     $y'_i = x_i \sin(-\theta^*) + y_i \cos(-\theta^*)$
Fit standard quadratic polynomial $y' = c_2 (x')^2 + c_1 x' + c_0$ to $(x'_i, y'_i)_{i=1}^N$ via least squares
Compute vertex in the rotated frame:
     $x'_v = -c_1 / (2 c_2)$
     $y'_v = c_0 - c_1^2 / (4 c_2)$
Derotate vertex back to original Cartesian frame:
     $x_v = x'_v \cos(\theta^*) - y'_v \sin(\theta^*)$
     $y_v = x'_v \sin(\theta^*) + y'_v \cos(\theta^*)$
Convert polynomial coefficients $[c_2, c_1, c_0]$ and rotation angle $\theta^*$ to algebraic coefficients $[A, B, C, D, E, F]$
return $(x_v, y_v)$ , $\theta^*$ , $[A, B, C, D, E, F]$ , $[c_2, c_1, c_0]$ , cost
Knowl 5 — Direct Least Squares Ellipse Initialization
model/method
bootEllipse generates an algebraic direct least-squares estimate of an ellipse from 2D coordinates $(x_i, y_i)_{i=1}^N$ . Based on the direct least-squares ellipse fitting formulation (Fitzgibbon-Pilu-Fisher and Halir-Flusser), it solves the constrained linear system: $\min_{\mathbf{a}} \|\mathbf{D} \mathbf{a}\|^2 \quad \text{subject to} \quad \mathbf{a}^T \mathbf{C} \mathbf{a} = 1$ where $\mathbf{a} = [A, B, C, D, E, F]^T$ , $\mathbf{D}$ is the $N \times 6$ design matrix with rows $[x_i^2, x_i y_i, y_i^2, x_i, y_i, 1]$ , and $\mathbf{C}$ is the $6 \times 6$ constraint matrix enforcing the ellipse condition $4AC - B^2 = 1$ .

The function outputs the estimated general quadratic coefficients $\mathbf{p}_A = [A, B, C, D, E, F]$ and the estimated centroid $(x_c, y_c)$ , which serve as initial parameters for non-linear optimization in fitConic.
Knowl 6 — Hyperbola Bootstrap Estimation
model/method
bootHyperbola generates an initial approximate parameter set for a hyperbola from 2D data points $(x_i, y_i)_{i=1}^N$ to bootstrap non-linear optimization in fitConic.

The algorithm minimizes a figure-of-merit cost function (fhyp, fhypopt) over a 3-parameter hyperbolic formulation and an initial rotation angle $\theta_0$ using numerical optimization (optim / optimize) up to a maximum iteration cap (maxiter, default $10000$ ). It returns:
- $\mathbf{p}_A$ : Six-element algebraic quadratic coefficient vector $[A, B, C, D, E, F]$ for the non-rotated hyperbola.
- $\mathbf{p}_{Ar}$ : Six-element algebraic quadratic coefficient vector for the rotated hyperbola.
- $\theta$ : The estimated rotation angle between $\mathbf{p}_A$ and $\mathbf{p}_{Ar}$ .
- fitdat: Optimization diagnostics and metadata returned by optim.
Knowl 7 — Synthetic Conic Section Data Generation
algorithm
The createConic procedure generates an array of $(x, y)$ coordinate pairs along a specified conic section over an evaluation grid $x \in \mathbb{R}^K$ , with optional random noise injection.
Input: Abscissa vector $X = (x_k)_{k=1}^K$ , parameters $\mathbf{p}$ (either $[A,B,C,D,E,F]$ or $[h,v,a,b,\theta]$ ), conic type $T \in \{\text{'e'}, \text{'h'}\}$ , noise generator $\text{ranFun}$ (optional), noise magnitude $\sigma \in \mathbb{R}$ (default 1), tolerance $\text{tol}$ (default $10^{-6}$ )
Output: $N \times 2$ array of $(x, y)$ points on the conic curve
if $\mathbf{p}$ has 5 geometric elements $[h,v,a,b,\theta]$ then
    Convert $\mathbf{p}$ to algebraic parameters $[A,B,C,D,E,F]$ using conic type $T$
end if
for each $x_k \in X$ do
    Solve $C y^2 + (B x_k + E) y + (A x_k^2 + D x_k + F) = 0$ for real roots $y$ :
     $\Delta_y = (B x_k + E)^2 - 4 C (A x_k^2 + D x_k + F)$
    if $\Delta_y \ge 0$ then
         $y^{(1)} = (-(B x_k + E) + \sqrt{\Delta_y}) / (2C)$
         $y^{(2)} = (-(B x_k + E) - \sqrt{\Delta_y}) / (2C)$
        Append $(x_k, y^{(1)})$ and $(x_k, y^{(2)})$ to output array
    end if
end for
if $\text{ranFun}$ is provided then
    for each generated point $(x_j, y_j)$ do
         $y_j \leftarrow y_j + \text{ranFun}(y_j) \times \sigma$
    end for
end if
return point array
Knowl 8 — Integer Observation Weighting for Conic Fitting
model/method
doWeights incorporates observation weights into 2D point data $\mathbf{XY} \in \mathbb{R}^{N \times 2}$ for weighted conic fitting.

Given an input weight vector $\mathbf{w} \in \mathbb{R}^N$ :
1. Negative weights are clamped to zero ( $w_j \leftarrow 0$ for $w_j < 0$ ).
2. Floating-point values are rounded to the nearest non-negative integer ( $w_j \leftarrow \lfloor w_j + 0.5 \rfloor$ ).
3. Row $j$ of $\mathbf{XY}$ is duplicated exactly $w_j$ times in the resulting dataset. Data points assigned zero weight ( $w_j = 0$ ) are removed entirely.
Knowl 9 — Conic Jacobian and Orthogonal Projection Residuals
model/method
During Levenberg-Marquardt fitting in fitConic, geometric error and derivatives are computed across iterations:
- Residual and Projection Calculation (Residuals.ellipse, Residuals.hyperbola): For current geometric parameters $\mathbf{p}_G = [h, v, a, b, \theta]$ and input points $\mathbf{XY} \in \mathbb{R}^{N \times 2}$ , the function computes the projected coordinates $\mathbf{XY}_{\text{proj}} \in \mathbb{R}^{N \times 2}$ on the conic curve and the root-sum-square (RSS) figure of merit: $\text{RSS} = \|\mathbf{XY} - \mathbf{XY}_{\text{proj}}\|_2 = \sqrt{\sum_{i=1}^N \left((x_i - x_{i,\text{proj}})^2 + (y_i - y_{i,\text{proj}})^2\right)}$
- Jacobian Matrix Construction (JmatrixLMA): Given algebraic parameters $\mathbf{p}_A = [A, B, C, D, E, F]$ and projected points $\mathbf{XY}_{\text{proj}}$ , it evaluates the $N \times 6$ Jacobian matrix $\mathbf{J}$ with columns corresponding to the general quadratic terms $[x^2, xy, y^2, x, y, 1]$ for each point, together with the residual norm vector $\mathbf{r} = \|\mathbf{XY} - \mathbf{XY}_{\text{proj}}\|$ .

Coverage note — No substantial contributed material was omitted; all core algorithms, fitting workflows, algebraic conversions, bootstrapping estimators, and geometric transformations in the package documentation have been fully covered.

References

1.https://www.mathworks.com/matlabcentral/answers/80541 for the RANSAC-style search to fit rotated parabolas. https://math.stackexchange.com/questions/426150 for detailed ellipse parametric equations. https://math.stackexchange.com/questions/2800817 for "focus/directrix/eccentricity" information https://people.cas.uab.edu/~mosya/cl/ and the folks referred to there, for fitConicLMA . https://en.wikipedia.org/wiki/Ellipse for several parameter conversion formulas
2.A. W. Fitzgibbon, M. Pilu, R. B. Fisher, "Direct Least Squares Fitting of Ellipses", IEEE Trans. PAMI, Vol. 21, pages 476-480 (1999)
3.Halir R, Flusser J (1998) Proceedings of the 6th International Conference in Central Europe on Computer Graphics and Visualization, Numerically stable direct least squares fitting of ellipses (WSCG, Plzen, Czech Republic), pp 125132.

Citation

MLA

Fitzgibbon, A., et al. “Direct Least Square Fitting of Ellipses”. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 21, no. 5, 1999, pp. 476–80, https://doi.org/10.1109/34.765658.

APA

Fitzgibbon, A., Pilu, M., & Fisher, R. B. (1999). Direct least square fitting of ellipses. IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(5), 476–480. https://doi.org/10.1109/34.765658

Chicago

Fitzgibbon, A., M. Pilu, and R. B. Fisher. 1999. “Direct Least Square Fitting of Ellipses”. IEEE Transactions on Pattern Analysis and Machine Intelligence 21 (5): 476–80. https://doi.org/10.1109/34.765658.

Harvard

Fitzgibbon, A., Pilu, M. and Fisher, R.B. (1999) “Direct least square fitting of ellipses”, IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(5), pp. 476–480. Available at: https://doi.org/10.1109/34.765658.

Vancouver

1. Fitzgibbon A, Pilu M, Fisher RB (1999) Direct least square fitting of ellipses. IEEE Transactions on Pattern Analysis and Machine Intelligence 21:476–480

BibTeX

@article{Fitzgibbon_1999, title={Direct least square fitting of ellipses}, volume={21}, ISSN={0162-8828}, url={http://dx.doi.org/10.1109/34.765658}, DOI={10.1109/34.765658}, number={5}, journal={IEEE Transactions on Pattern Analysis and Machine Intelligence}, publisher={Institute of Electrical and Electronics Engineers (IEEE)}, author={Fitzgibbon, A. and Pilu, M. and Fisher, R.B.}, year={1999}, month=May, pages={476–480} }

Metadata:Crossref

Access the Paper

This paper is available from its original source. Click below to access the PDF.

Open PDF

Direct Least Square Fitting of Ellipses

Abstract

Topics

Keywords

Table of Contents

Knowls

Knowl 1 — Conic Section Fitting via Levenberg-Marquardt Optimization

Knowl 2 — Algebraic and Geometric Parameterizations of Conic Sections

Knowl 3 — Derotation Angle Formula for Quadratic Conics

Knowl 4 — Rotated Parabola Fitting Algorithm

Knowl 5 — Direct Least Squares Ellipse Initialization

Knowl 6 — Hyperbola Bootstrap Estimation

Knowl 7 — Synthetic Conic Section Data Generation

Knowl 8 — Integer Observation Weighting for Conic Fitting

Knowl 9 — Conic Jacobian and Orthogonal Projection Residuals

References

Citation

BibTeX

Access the Paper