14th IAPR International Conference on Machine Vision Applications (MVA)

May 18-22, 2015. Miraikan, Tokyo, Japan

Grouped Outlier Removal for Robust Ellipse Fitting

Mang Shao Yoshihisa Ijiri and Kosuke Hattori
Imperial College London OMRON Corporation
[email protected] [email protected]
kosuke [email protected]

Abstract
This paper presents a novel outlier removal method
which is capable of fitting ellipse in real-time under
high outlier rate, based on the phenomenon that out-
liers generated by ellipse edge point detector are likely
to appear as groups due to real-world nuisances, such
as under partial occlusion or illumination change. To Figure 1. A challenging case for proximity-based
confront the grouped outliers while maintaining the fit- outlier removal method. Type (a) outliers cannot
ting efficiency, we introduce a proximity-based ‘split be filtered by simple proximity (e.g. k-Nearest
and merge’ approach to cluster the edge points into sub- Neighbour) check. Type (b) is even more difficult
sets, followed by a breath-first outlier removal process. since they are connected with the inlier contour.
The experiment shows that our algorithm achieves high
performance under a wide range of inlier ratio and
noise level with various types of realistic nuisances.
type (b). Another work [8] presents an accurate, non-
1 Introduction iterative method based on the geometric distance be-
tween a data point and an ellipse. However it does not
take the type (b) outlier into account as well.
Ellipse fitting is one of the fundamental problems in In this work, we demonstrate how edge points can
computer vision and robotic tasks. It is required as be effectively grouped into short contours to reduce the
preprocessing modules in many higher level applica- computational cost, and further show how to eliminate
tions such as textureless object recognition and shape the outlier contours in a breath-first manner. Our con-
alignment. In this paper, we propose a real-time out- tribution is four-fold. First we introduce a split-and-
lier removal solution for ellipse fitting which efficiently merge trick to cluster data points into subsets which
deals with outliers especially when they are contami- are likely to contain either only inliers or only outliers.
nated consistently. Second we propose a breath-first strategy for search-
Amongst ellipse fitting algorithms presented in the ing the outlier contours through the combination of
last few decades, efforts have been made to approach subsets. Third, we speed up the searching process by
the KCR lower bound [1], i.e. the theoretical accu- using the smallest generalised eigenvalue (which is a
racy bound for ellipse fitting without outliers. In these by-product of algebraic fitting algorithm) to approxi-
works, the presence of noises is considered from ob- mate the point-to-curve projection, and then use the
servational error only and thus the point-to-curve pro- algebraic fitting solution as initial guess for geometric
jection error is mostly assumed under Gaussian dis- fitting algorithm to alleviate measurement noises. Fi-
tribution. While in the real-world scenario, image or nally, we proposed a synthetic dataset for evaluation.
extracted edgess is likely to be contaminated in a very
biased way such as under partial occlusion, specular
highlight, deformation or shading. In such cases, ro- 2 Preliminaries and Related Works
bust fitting algorithm like random sample consensus
(RANSAC) [4] is generally applied to eliminate out- Given a point set, the objective of ellipse fitting
liers. is to find a geometric parameter set that minimises
When the inlier rate  is low, RANSAC soon be- the sum of inlier-points-to-ellipse-curve projection dis-
comes infeasible for many applications since the possi- tance. In this section, we briefly introduce two broad
bility of finding at least one correct ellipse model after algorithm categories: algebraic (Least-Square-based)
K iterations p = 1 − (1 − 5 )K decreases drastically. and geometric (Maximum-Likelihood-based) method.
Also, the model candidates generated from minimal
sample sets can be greatly affected by the noise, which 2.1 Algebraic Fitting
leads to suboptimal solutions.
A recent work [10] proposed a proximity-based out- Any ellipse can be represented by a second order
lier detection algorithm to effectively remove the iso- polynomial F (u, x) = u · x = ax2 + bxy + cy 2 + dx +
lated outliers and outlier clusters by constructing a ey + f = 0, subject to b2 < 2ac. Our goal here is
proximity graph. However, the adjacency matrix is to estimate the parameter set u = [a, b, c, d, e, f ] from
expensive to compute if data point set is large and a given point set xi = [xi , yi ] such that the sum of
some parameters need to be tuned carefully in or- 
N
algebraic distance |F (u, xi )|2 is minimised.
der to achieve a good clustering result. Also, the i=1
proximity-based method fails if the outliers are con- This problem is generally tackled with linear least
nected smoothly with the inliers, as shown in Figure 1, square solvers as in several seminal approaches [9, 5, 6].
978-4-901122-14-6 ©2015 MVA Organization 138
In this work, since we focus on outlier elimination, and a adequate connection radius is also needed to be
Taubin’s method [9] is chosen as our algebraic ellipse carefully chosen. The edge point detector we used in
fitting algorithm. From our empirical experiments, we this work naturally provides the connectivity between
believe that Taubin’s method is still one of the most ac- the points, which greatly simplifies the clustering pro-
curate and robust methods giving its efficiency, despite cess.
it has been proposed for decades. Since this method
is designed for general conic fitting, it might return Algorithm 1: Proximity-based Point Clustering
conics other than ellipse as well (e.g. hyperbola curve)
Input: Cf ull = {xi }N i=1
and hence leads to incorrect fitting result.
Initialisation: t, τ, D, C1 = {};
for i ← 1 to N do
2.2 Geometric Fitting Ck = Ck ∪ x i ;
di = |xi − xi+1 |2 ; /* i + 1 = 1 if i = N */
Geometric parameters of ellipse is of 4 elements: cen- if di > t ∗ median{di } then
i
ter position, length of major axis, minor axis and angle k = k + 1; Ck = {};
of tilt. Unlike algebraic methods, the linear condition end
has no longer held for solving the least square prob- end
lem with geometric parameters. Instead, iterative op- Delete any set Ck that has cardinality |Ck | < τ
for k ← 1 to K do
timisation methods are used to find the local optimal
if |Ck | > 2N
D
then
solution given an initial guess. |C |
Uniformly split Ck into  rk  sets;
Geometric fitting, i.e. Maximum-Likelihood-based Replace Ck by these sets;
method (ML) is generally regarded as one of the most end
precise fitting algorithms. They do not suffer from else if |Ck | + |Ck+1 | < N then
D
scale indeterminacy as in the algebraic methods and Ck = Ck ∪ Ck+1 ; /* Delete set Ck+1 */
are likely to achieve the local optimal solution. A de- end
tailed analysis is given in the work [6]. end
The main drawback of geometric fitting algorithms return {Ck }K k=1
is their high computational cost due to iterative op-
timisation process. With existence of outliers or high
measurement noise, they also tend to fail to converge. The point set is ordered by the edge point detector
The initial guess is crucial as well to achieve the global such that xn and xn+1 are next to each other. A point
optimal solution. set is split at a point that Euclidean distance to its
In this work, we apply ML as the final refinement neighbour points di greater than t ∗ median{di }, where
i
to alleviate the measurement noise while keeping the t is the distance ratio threshold. If the measurement
pipeline fast by using a near-optimal algebraic solution noise is too high, t should be set to a larger number
as its initial guess. to prevent over-segmentation. In this step, subset is
discarded as isolated outliers if it contains less than τ
2.3 Edge Point Detector points.
To deal with type (b) outliers in Figure 1, each sub-
In this work, the edge points are collected circularly set is uniformly split if they contain more than a cer-
around a user-estimated center point by finding the tain number of points. This number is determined by
maximal gradient change along each circle radius such a pre-defined expected subset cardinality D. Larger D
that the order of the point sequence is naturally known. allows a finer segmentation of edge points, but induces
The center point does not necessarily need to be accu- higher risk of getting stuck in local optimum after the
rate. This setting is common in many industrial ap- later searching stage, and also sacrifices the processing
plications like shape alignment, as it provides several speed.
merits over the ordinary edge detector.
In the last step, neighbouring subsets are merged if
their sum of cardinality is sufficiently small. In the
3 Proposed Approach end, the whole point set should be split into similar
sizes and each subset is likely to contain only inliers or
The proposed method consists of three stages: outliers. The pseudo code of this section is shown in
• Clustering data points based on proximity so that the Algorithm 1.
each subset is likely to contain only inliers or out-
liers edge points. 3.2 Breath-First Searching
• Searching through combinations of subsets to min- As shown in Figure 2, the fitting trials are per-
imise algebraic fitting distance until convergence. formed between the combinations of subsets only.
The total number of possible combinations {C}D n=1 is
• Refining the algebraic solution by the geometric

D
D!
fitting algorithm.
n!(D−n)! , which could be still large in number.
n=0
3.1 Proximity-based Point Clustering For each test, Taubin’s fitting method [9] is used
to sidestep the expensive geometric projection. We
Generally a partial or full adjacent matrix needs to further speed up the process by approximating alge-
be calculated to determine the connections between braic distance with the smallest generalised eigenvalue
edge points. This process roughly has computational λ from solving the linear least square system.
complexity O(N 2 ) respect to the number of data points The algebraic distance between an ellipse u and a
139
 Figure 2. This figure shows an example that the value of the smallest generalised eigenvalue λ decreasing
drastically after excluding the outlier subsets. The point set contains two ‘outlier‘ subsets, (a), (b) and (c)
show three cases during the iteration, while (a) does not exclude any outlier subset, (b) excludes one of them
and (c) excludes all.

point set C is: Algorithm 2: Preemptive searching

 Input: {Ck }K
F (u, C) = |u · Z(x)|2 , Z(x) = [x2 , xy, y 2 , x, y, 1]. k=1
Initialisation: S, σ, η, {Os = ∅}Ss=1 ;
x∈C for m ← 1 to  K do
Ctest = k∈K,k=m {Ck } ;
We define our energy function to approximate the al- [λm,1 , um,1 ] = T aubinF itting(Ctest );
gebraic distance by the minimal generalised eigenvalue end
solved as follows: Os = {m}, Es = λm for sth least λm ;
for l ← K − 2 to 1 do
E(C) = min(λ), subject to Mv = λNv for s ← 1 to |{Os }| do
for m ← 1 to  K, m ∈ Os do
N
1  Ctest = k∈K,k=m {Ck } ;
where M = Z(xT − x̄T )Z(x − x̄) [λm,s , um,s ] = T aubinF itting(Ctest );
N n=1
end
⎛ ⎞ end
x2n x n yn 0 f0 xn 0 0
⎜x n y n x2n + yn2 x n yn f 0 yn f0 xn 0⎟ if (E1 > σ ∗ min {λm,s })&(E1 < η) then
4
|C| ⎜
⎜ 0 x n yn yn2 0 f 0 yn 0⎟
⎟ return uarg min{λm,s }
N = |C| n=1 ⎜ f x ⎟ (m,s)
⎜ 0 n f 0 yn 0 f02 0 0⎟ end
⎝ 0 f0 xn f 0 yn 0 f02 0⎠ else
0 0 0 0 0 0 S ← l if S > l;
Here we propose a searching scheme that efficiently for s ← 1 to |{Os }| do
replace {Os } by {Os ∪ {m}} for S least λm,s ;
search and exclude the outlier subsets in a breath-first end
manner. To minimise the energy while maximising update Es ;
the number of possible inliers, the combination of sets end
are tested by excluding one subset in each iteration. end
The searching process converges if excluding any sub-
set does not reduce the energy up to certain rate σ. To
prevent local optimal solution, number of S candidate
sets with lowest energy are stored at each iteration, to the optimal solution. Thus the efficiency is main-
and each will be appended with another S candidates tained despite the algorithm itself is expensive because
in the next iteration. In case there are much more it converges within few steps in most of cases.
outliers than average subset size so that removing any
subset does not reduce the energy, we introduce an er- 4 Evaluation
ror threshold η as another stop condition. The pseudo
code of this section is provided in Algorithm 2. Our method is evaluated by both synthetic1 and re-
alistic datasets which consists of ellipse edge points
3.3 Refinement with different types of contamination. For synthetic
data, we have rendered 3 videos (2 with occlusions and
Experimentally we found a refinement step to be cru- 1 with background clutter) with 3D computer graphic
cial to reducing the measurement noise from the edge software at 960 × 540 resolution. Each video contains
point detector. Given the inlier point set, we convert 100 frames and each frame consists of 100 extracted
the algebraic parameters u to geometric parameters H ellipse edge points. The ground truth ellipse centers
and apply a nonlinear optimiser as below to further are fixed to the image center. To simulate the mea-
minimise the point-to-curve (orthogonal to the ellipse surement noise from camera, we augmented Gaussian
tangent) projection error. The detail for calculating noise (σ = 3 pixels) to the data points. We use the
projection distance can be found in [3]. following parameter setting for our method in all ex-
periments: t = 2, τ = 5, D = 12, S = D/2, η = 5 and
H ∗ = arg min ||proj(H, C)||2 σ = 1.05. The preliminary investigation suggests this
H
parameter set adapts a wide range of inlier rate and
Since the initial geometric parameters are already es-
timated from the previous stages, they are close enough 1 http://bit.ly/1Dvs0id

140
 Clutter 1
Err < 5 pix Clutter 1 Occ. 1 Occ. 2 Taubin
Taubin 45.97% 50.97% 2.78% Taubin + RANSAC
+ RANSAC 66.17% 21.78% 23.29% Our result
Our result 93.80% 85.22% 92.27% + Edge points
Clutter 1
Occlusion 1

Occlusion 2
Occlusion 1

Synthetic dataset (quantitative)

Our result
Occlusion 2
Edge points

*Constrained to circle
Synthetic dataset (qualitative) Realistic dataset (qualitative)

Figure 3. Evaluation result of 3 methods on synthetic and realistic datasets (best view in color). The boxplot
reflects the statistic of center errors (capped at 20 pixel) from 100 runs of each method on every video frames.
The table on the top left corner shows the success rate (center error less than 5 pixels) of the methods on
each video. Note that for synthetic dataset, Gaussian noise is augmented to all data points in each run to
simulate the noise from realistic optical sensors.

outlier type. which is sufficient to achieve real-time processing speed

We compare our method with Taubin’s method for many higher level applications. Finally, we believe
with/without RANSAC and summarise the evaluation that our approach extends naturally to general conic
result in Fig 3. In overall, our approach achieves supe- and ellipsoid fitting.
rior accuracy than the baselines.
For Taubin’s method itself, almost all fitting esti- References
mations are distorted by outliers due to the nature of
least square solvers, such as in the dataset ‘Occlusion 2’ [1] Nikolai Chernov and Claire Lesort. Statistical effi-
shown in Fig 3. However, RANSAC does not boost the ciency of curve fitting algorithms. Computational statis-
accuracy significantly as expected as well. The main tics & data analysis, 47(4):713–728, 2004.
reason is twofold. Firstly the parameters, e.g. inlier [2] Ondřej Chum, Jiřı́ Matas, and Josef Kittler. Locally
ratio and distance threshold, are employed empirically optimized ransac. In Pattern Recognition, pages 236–
as stop criteria. Such settings are not optimal to all 243. Springer, 2003.
situations therefore leading to an overall poor perfor- [3] David Eberly. Distance from a point to an ellipse, an
mance. Another is due to the high measurement noise ellipsoid, or a hyperellipsoid. Technical report, Tech-
within all data points. Since RANSAC picks hypothe- nical report, Geometric Tools, LLC, 2011. 5.1. 3, 1998.
ses from random minimal samples, the impact from [4] Martin A Fischler and Robert C Bolles. Random sam-
noises is drastically amplified comparing to estimating ple consensus: a paradigm for model fitting with appli-
from all inlier samples, which also explains the mo- cations to image analysis and automated cartography.
tivation of our method. There are several extended Communications of the ACM, 24(6):381–395, 1981.
RANSAC-like methods have been proposed to deal [5] Andrew Fitzgibbon, Maurizio Pilu, and Robert B Fisher.
with such a problem, they either are expensive (e.g. Direct least square fitting of ellipses. Pattern Analy-
pre-emptive RANSAC[7]) for ellipse fitting problem or sis and Machine Intelligence, IEEE Transactions on,
suffering from poor initial estimations generated from 21(5):476–480, 1999.
minimal samples(e.g. locally optimised RANSAC[2]). [6] Kenichi Kanatani and Prasanna Rangarajan. Hyper-
Our method achieves less than 20ms runtime on accurate ellipse fitting without iterations. In VISAPP
single core CPU with up to 360 data points, which (2), pages 5–12, 2010.
meets the time requirement as a pre-processing mod- [7] David Nistér. Preemptive ransac for live structure and
ule for many higher level real-time applications. With motion estimation. Machine Vision and Applications,
code optimisation and parallel processing, the whole 16(5):321–329, 2005.
pipeline can be further speeded up. [8] Dilip K Prasad, Maylor KH Leung, and Chai Quek.
Ellifit: An unconstrained, non-iterative, least squares
5 Conclusion based geometric ellipse fitting method. Pattern Recog-
nition, 46(5):1449–1465, 2013.
This paper has presented a new approach for effi- [9] Gabriel Taubin. Estimation of planar curves, surfaces,
cient ellipses fitting under high outlier rate. We have and nonplanar space curves defined by implicit equa-
demonstrated how data points can be clustered based tions with applications to edge and range image seg-
on their proximity and how outliers can be filtered in mentation. IEEE Transactions on Pattern Analysis
a ‘preemptive’ manner. Our method has shown to be and Machine Intelligence, 13(11):1115–1138, 1991.
especially effective on the ‘grouped’ outliers due to en- [10] Jieqi Yu, Haipeng Zheng, Sanjeev R Kulkarni, and
vironmental nuisances such as partial occlusion, spec- H Vincent Poor. Two-stage outlier elimination for ro-
ular highlight, deformation or shading. The runtime bust curve and surface fitting. EURASIP Journal on
of the method is below 20ms for up to 360 data points, Advances in Signal Processing, 2010:4, 2010.

141