Expert Systems With Applications 222 (2023) 119853

Contents lists available at ScienceDirect

Expert Systems With Applications

journal homepage: www.elsevier.com/locate/eswa

A new efficient method for solving the multiple ellipse detection problem
Rudolf Scitovski ∗, Kristian Sabo, Patrick Nikić, Snježana Majstorović
Department of Mathematics, University of Osijek, Trg Ljudevita Gaja 6, HR - 31000 Osijek, Croatia

ARTICLE INFO ABSTRACT

Keywords: In this paper, we consider the multiple ellipse detection problem based on data points coming from a number
Multiple ellipse detection problem of ellipses in the plane not known in advance. In so doing, data points are usually contaminated with some
Mahalanobis k-means algorithm noisy errors. In this paper, the multiple ellipse detection problem is solved as a center-based problem from
Global optimal partition
cluster analysis. Therefore, an ellipse is considered a Mahalanobis circle. In this way, we easily determine a
DIRECT
distance from a point to the ellipse and also an ellipse as the cluster center. In the case when the number
of ellipses is known in advance, an optimal partition is searched for on the basis of the 𝑘-means algorithm
that is modified for this case. Hence, a good initial approximation for M-circle-centers is searched for as unit
circles with the application of a few iterations of the well-known DIRECT algorithm for global optimization.
In the case when the number of ellipses is not known in advance, optimal partitions with 1, 2, … clusters for
the case when cluster-centers are ellipses are determined by using an incremental algorithm. Among them,
the partition with the most appropriate number of clusters is selected. For that purpose, a new Geometrical
Objects-index (GO-index) is defined. Numerous test-examples point to high efficiency of the proposed method.
Many algorithms can be found in the literature that recognize ellipses with clear edges well, but that do
not recognize ellipses with unclear or noisy edges. On the other hand, our algorithm is specifically used for
recognition of ellipses with unclear or noisy edges.

1. Introduction processing in industrial applications (Akinlar & Topal, 2013; Prasad

et al., 2013), astronomical and geological shape segmentation (Wynn &
We consider the multiple ellipse detection problem on the basis Stewart, 2005), wireless sensor networks (Moshtaghi et al., 2011), etc.
of the set of data points coming from a number of ellipses in the Many studies can be found in the literature that deal with this issue
plane not known in advance. Data points are usually contaminated with in different ways. We will list some that are frequently mentioned and
some noisy errors. A new method proposed in this paper recognizes compare them with our method. The algorithm proposed in Akinlar
ellipses with unclear or noisy edges in a very short time. The method and Topal (2013) recognizes ellipses with clear edges in a very short
is based on center-based clustering, where an ellipse is interpreted time, but it recognizes neither ellipses with unclear or noisy edges nor
as a Mahalanobis circle. The solution is searched for by means of a ellipses with a semi-axes ratio greater than 4. The algorithm proposed
modification of the well-known 𝑘-means algorithm with a good initial in Fornaciari et al. (2014) relies on an innovative selection strategy
approximation obtained in advance by using the global optimization of arcs which are candidates to form ellipses and on the use of the
algorithm DIRECT (see e.g. Jones & Martins, 2020; Jones et al., 1993). Hough transform to estimate parameters in a decomposed space. As
A new index was constructed for the purpose of recognizing the par- shown in Scitovski and Marošević (2014), for solving such problems,
tition with the most appropriate number of clusters (MAPart). In a the method based on the Hough transformation is inferior in relation
large number of cases with data that do not contain outliers, the new to the incremental algorithm, and the incremental algorithm performed
index shows the property of monotonicity, i.e., with an increasing worse in relation to the statistically based method proposed by Grbić
number of clusters its value decreases monotonically to the MAPart. et al. (2016). The algorithm proposed by Grbić et al. (2016) combines
The use of the monotonicity of the new index significantly accelerates very successfully the well-known direct least squares method and the
the algorithm. RANSAC algorithm with a realistic statistical model of multiple ellipses
This problem occurs in numerous applications. Let us mention only in the plane. In similar examples, the new proposed method shows sig-
a few, e.g., medical diagnosis (Rueda et al., 2013), pattern recogni- nificant improvements. An unconstrained, non-iterative, least squares
tion and computer vision (Liu & Ramani, 2009), agriculture (Kashiha based geometric ellipse fitting method was proposed in Prasad et al.
et al., 2014), robotics, object detection, pupil tracking, and other image

∗ Corresponding author.
E-mail addresses: [email protected] (R. Scitovski), [email protected] (K. Sabo), [email protected] (P. Nikić), [email protected] (S. Majstorović).

https://doi.org/10.1016/j.eswa.2023.119853
Received 11 November 2022; Received in revised form 19 January 2023; Accepted 8 March 2023
Available online 11 March 2023
0957-4174/© 2023 Elsevier Ltd. All rights reserved.
R. Scitovski et al. Expert Systems With Applications 222 (2023) 119853

(2013), and fast ellipse recognition based on the three-point algorithm ellipses as Mahalanobis circle-centers. Section 4 describes a data set
with edge angle information was proposed in Kwon et al. (2016). A and refers to our publicly available software. This section also provides
fast and robust ellipse detection method is proposed in Dong et al. the results of algorithm testing. Finally, some conclusions are given in
(2018). The method adopts arcs selection, smart grouping and gradient Section 5.
information. In the paper by Liu et al. (2022), the method for ellipse
detection from images, based on a deep learning model that can intelli-
gently distinguish between useful edges and useless edges. In the paper 2. An ellipse detection problem
by Lu et al. (2020), the authors propose a valuable industry-oriented
ellipse detector by arc-support line segments, which simultaneously Let  = {𝑎𝑖 = (𝑥𝑖 , 𝑦𝑖 )𝑇 ∈ R2 ∶ 𝑖 = 1, … , 𝑚},  ⊂ 𝛥 ∶= [𝛼1 , 𝛽1 ]×[𝛼2 , 𝛽2 ],
reaches high detection accuracy and efficiency. A fast and robust ellipse be a set of data points coming from an ellipse in the plane not known
detection method based on sorted merging is proposed by Wang et al. in advance that should be reconstructed or detected.
(2014). This method first represents the edge bitmap approximately by
a set of line segments, and then gradually merges the line segments
into elliptical arcs and ellipses. A method for automatic and robust 2.1. Searching for an ellipse in classical form
detection of non-overlapping ellipses from 2D points is proposed in
the paper by Maalek and Lichti (2021). This method includes three In general, an ellipse in the plane can be defined as a set of points
new developments: (i) selecting initial subsamples using a bucketing {(𝑥, 𝑦) ∈ R2 ∶ 𝑃 (𝑥, 𝑦) = 0}, where
strategy based on the polar angle of the points; (ii) detecting inlier
points by reducing the robust ellipse fitting to a robust circle fitting 𝑃 (𝑥, 𝑦) = 𝐴𝑥2 + 2𝐵𝑥𝑦 + 𝐶𝑦2 + 𝐷𝑥 + 𝐸𝑦 + 𝐹 , 𝐴𝐶 − 𝐵 2 > 0. (1)
problem; and (iii) choosing the best inlier set among all subsamples √
2
Dividing the equation 𝑃 (𝑥, 𝑦) = 0 by 𝐴𝐶 − 𝐵 , we obtain the equation
using adaptive, systematic, and generalizable selection criteria. The
paper by Liu et al. (2019) presents a fast and robust ellipse detector.
The detector first extracts segments using an edge predictor based on 𝑎𝑥2 + 2𝑏𝑥𝑦 + 𝑐𝑦2 + 𝑑𝑥 + 𝑒𝑦 + 𝑓 = 0, (2)
curvature analysis. Then, line segments are generated based on the
length condition other than least squares approximation. After that, where 𝑎𝑐 − 𝑏2 = 1. Now, for a given set  we will search for parameters
potential ellipses are detected based on edge curvature and convexity. 𝑎, 𝑏, 𝑐, 𝑑, 𝑒, 𝑓 ∈ R, with 𝑎𝑐 − 𝑏2 = 1, by solving a linear least squares
In addition, a re-find contours detection method is introduced to im- problem
prove the accuracy by searching edge points in the missing part of the
∑
𝑚
ellipse. In Jia et al. (2017), a new geometry constraint is introduced to argmin (𝑎𝑥2𝑖 + 2𝑏𝑥𝑖 𝑦𝑖 + 𝑐𝑦2𝑖 + 𝑑𝑥𝑖 + 𝑒𝑦𝑖 + 𝑓 )2 . (3)
prune the undesired candidates and to pick out elliptical ones. Based 𝑎,𝑏,𝑐,𝑑,𝑒,𝑓 ∈R
𝑖=1
𝑎𝑐−𝑏2 =1
on the projective invariant, it is able to reflect the intrinsic geometry
of a planar curve. Some other approaches to solving this problem can be found in Fitzgib-
Unlike the aforementioned algorithms that recognize well ellipses bon et al. (1999).
with clear edges, but do not recognize ellipses with unclear or noisy
edges, our algorithm is specifically used for recognition of ellipses with 2.2. Searching for an ellipse in trigonometric form
unclear or noisy edges. Furthermore, since our method is based on
center-based clustering and an incremental algorithm, where an ellipse
is interpreted as a Mahalanobis circle, it is extremely important to By knowing the values of parameters 𝑎, 𝑏, 𝑐, 𝑑, 𝑒, 𝑓 , ellipse 𝐸 can also
recognize the MAPart. For this purpose, a specialized index (a GO- be defined as a set of points (𝑥, 𝑦) ∈ R2 satisfying the equation (see
index) is proposed, for which numerous numerical experiments show e.g. Gander et al., 1994; Marošević & Scitovski, 2015; Prasad et al.,
that it recognizes the MAPartvery well. With appropriate modification, 2013):
this index can also be used for recognition of other geometrical objects, [(𝑥 − 𝑝) cos 𝜗 + (𝑦 − 𝑞) sin 𝜗]2 [(𝑥 − 𝑝) sin 𝜗 − (𝑦 − 𝑞) cos 𝜗]2
such as circles, lines, etc. (see e.g. Sabo et al., 2020; Scitovski & Sabo, + = 1, (4)
𝜉2 𝜂2
2019).
The main challenge and issue in this study is to detect multiple inter- where 𝑆 = (𝑝, 𝑞)𝑇 is the center of the ellipse, 𝜉, 𝜂 > 0 are the lengths
secting ellipses in the plane based on the set of data points coming from of semi-axes, and 𝜗 is the angle between semi-axes 𝜉 and a positive
these ellipses. That problem was successfully solved by interpreting an direction of the abscissa. Therefore, an ellipse will be denoted by
ellipse as the Mahalanobis circle. The multiple circle detection problem 𝐸(𝑆, 𝜉, 𝜂, 𝜗). [ ]
is much simpler (see Scitovski & Sabo, 2019). cos 𝜗 − sin 𝜗
Let 𝜗 ∈ [0, 𝜋], and let 𝑈 (𝜗) = be a rotation matrix.
The method proposed in this paper can be applied to multiple ellipse sin 𝜗 cos 𝜗
detection in real-world images. It is thereby necessary to do image Then the ellipse (4) can be written in the parametric form:
preprocessing first, which is based on some edge detection method [ ] [ ]
𝑥(𝑡) 𝜉 cos 𝑡
(see e.g. Mathematica computation system (Wolfram, 2020) or the = 𝑆 + 𝑈 (𝜗) , 𝑡 ∈ [0, 2𝜋]. (5)
𝑦(𝑡) 𝜂 sin 𝑡
Canny filter from OpenCV Bradski, 2000). Our method has proven
to be effective on numerous real-world images. Unfortunately, in the Now, for a given set of data points  we will search for parameters
case when the preprocessed image contains, in addition to ellipses, a 𝑝, 𝑞, 𝜉, 𝜂, 𝜗 by solving the following nonlinear least squares problem:
large number of non-ellipsoidal geometrical objects, the method shows ∑
𝑚
certain weaknesses compared to other previously described methods. argmin D(𝐸((𝑝, 𝑞), 𝜉, 𝜂, 𝜗), 𝑎𝑖 ), (6)
(𝑝,𝑞)∈𝛥
The paper is organized as follows. An ellipse detection problem is 𝜉,𝜂∈[0,𝑅],𝜗∈[0,𝜋]
𝑖=1
considered in Section 2. In particular, the case of an ellipse defined as
a Mahalanobis circle is observed in Section 2.3. The multiple ellipse where 𝑅 = 12 min{𝛼2 − 𝛼1 , 𝛽2 − 𝛽1 }. The distance D(𝐸, 𝑎𝑖 ) between the
detection problem is considered in Section 3. In particular, the case ellipse 𝐸 and the point 𝑎𝑖 ∈  can be determined numerically (see
when the number of ellipses is known in advance is observed in e.g. Uteshev & Goncharova, 2018), and the solution to GOP (6) can be
Section 3.1, and the case when the number of ellipses is not known in found by using the global optimization algorithm DIRECT (see Grbić
advance is observed in Section 3.2. In the latter case, a new, efficient et al., 2013; Jones & Martins, 2020; Jones et al., 1993). However, the
GO-index is proposed which enables recognition of the MAPart and efficiency of this method would be very low.

2
R. Scitovski et al. Expert Systems With Applications 222 (2023) 119853

2.3. Searching for an ellipse as a Mahalanobis circle where 𝑀2 is the set of positive definite symmetric matrices of order
two.
The ellipse 𝐸 given by (4) can be written in the following form: Solutions to problem (17) can be found by using some local opti-
mization method like Nelder–Meade or Quasi-Newton (see e.g. Dennis
𝐸 = {𝑢 ∈ R2 ∶ (𝑆 − 𝑢)𝑇 𝛴 −1 (𝑆 − 𝑢) = 1}, (7)
[ 2 ] & Schnabel, 1996; Nocedal & Wright, 2006) since we are able to de-
𝜉 0 √ termine a very favorable initial approximation M0 (𝑆 0 , 𝑟0 , 𝛴 0 ). Namely,
where 𝛴 = 𝑈 (𝜗) 𝑈 𝑇 (𝜗). If we choose 𝑟 ∶= 𝜉𝜂 as a geometric
0 𝜂2 for 𝑆 0 we can choose a centroid (mean) of the set . Furthermore,
√
mean of the semi-axes 𝜉 and 𝜂, i.e., 𝑟2 = det 𝛴, then Eq. (7) can be for the matrix 𝛴 0 we can choose a corresponding covariance matrix
∑
written as: 𝛴 0 = 𝑚1 𝑚 0 𝑖 0 𝑖 𝑇 0
𝑖=1 (𝑆 − 𝑎 )(𝑆 − 𝑎 ) , and 𝑟 is determined from
√
det 𝛴(𝑆 − 𝑢)𝑇 𝛴 −1 (𝑆 − 𝑢) = 𝑟2 . (8) ∑
𝑚
1
(𝑟0 )2 = 𝑚
‖𝑆 0 − 𝑎𝑖 ‖2 , (18)
This motivates us to give the following definition. 𝑖=1

because
Definition 1. Mahalanobis circle (M-circle) M(𝑆, 𝑟; 𝛴) with the center
∑
𝑚
( 0 )2 ∑
𝑚
( 0 )2
𝑆 ∈ R2 , the radius 𝑟, and the covariance matrix 𝛴 is the set of points ‖𝑆 − 𝑎𝑖 ‖2 − 𝑟2 ≥ ‖𝑆 − 𝑎𝑖 ‖2 − (𝑟0 )2 , for all 𝑟 ∈ R. (19)
𝑖=1 𝑖=1
M(𝑆, 𝑟; 𝛴) = {𝑢 ∈ R2 ∶ 𝑑𝑀 (𝑆, 𝑢; 𝛴) = 𝑟2 }, (9)
where the distance-like function 𝑑𝑀 ∶ R2 × R2 → R+ given by 3. Multiple ellipse detection problem
√
𝑑𝑀 (𝑆, 𝑢; 𝛴) = det 𝛴(𝑆 − 𝑢)𝑇 𝛴 −1 (𝑆 − 𝑢) = ‖𝑢 − 𝑣‖2𝛴 (10) The multiple ellipse detection problem will be considered as a
special center-based clustering problem, in the case when the number
is the Mahalanobis distance-like function.
of ellipses is known in advance and in the case when the number of
Especially, for 𝑟 = 1, set (9) is called a normalized Mahalanobis circle. ellipses is not known in advance.
In this case the product of semi-axes of the corresponding ellipse is It is well known that determining the most appropriate number of
equal to 1 and the area of this ellipse is equal to 𝜋. clusters is the most complex problem in cluster analysis.
In the sequel, we present a lemma that shows the relationship
between the ellipse 𝐸 and the M-circle M. Moreover, it gives an explicit 3.1. The number of ellipses is known in advance
formula for the transformation of the ellipse 𝐸 into an M-circle M and
vice versa. Let  be the set of data points coming from 𝑘 ≥ 1 ellipse in the plane
𝐸𝑗 (𝑆𝑗 , 𝜉𝑗 , 𝜂𝑗 , 𝜗𝑗 ), 𝑗 = 1, … , 𝑘, that should be reconstructed or detected.
Lemma 1. Ellipse 𝐸(𝑆, 𝜉, 𝜂, 𝜗) given√ by (4) or (5) can be considered an M- Since we consider the multiple ellipse detection problem as a special
circle M(𝑆, 𝑟; 𝛴), with radius 𝑟 = 𝜉𝜂. Conversely, an M-circle M(𝑆, 𝑟; 𝛴) center-based clustering problem, our problem is:
( 2 )
is an ellipse 𝐸(𝑆, 𝜉, 𝜂, 𝜗), where diag(𝜉 2 , 𝜂 2 ) = 𝑈 𝑇 (𝜗) √ 𝑟 𝛴 𝑈 (𝜗), 𝛴 =
[ ] [ ] det 𝛴 ∑
𝑚
𝜎11 𝜎12 cos 𝜗 − sin 𝜗 argmin 𝐹 (𝑆, 𝜉, 𝜂, 𝜗), 𝐹 (𝑆, 𝜉, 𝜂, 𝜗) = min D(𝐸𝑗 (𝑆𝑗 , 𝜉𝑗 , 𝜂𝑗 , 𝜗𝑗 ), 𝑎𝑖 ).
, 𝑈 (𝜗) = and 𝑆∈𝛥𝑘 1≤𝑗≤𝑘
𝜎12 𝜎22 sin 𝜗 cos 𝜗 𝜉,𝜂∈[0,𝑅]𝑘 ,𝜗∈[0,𝜋]𝑘
𝑖=1

⎧[0, 𝜋∕4), 𝜎12 ≥ 0 & 𝜎11 > 𝜎22

(20)
⎪
1 2𝜎12 ⎪[𝜋∕4, 𝜋∕2), 𝜎12 > 0 & 𝜎11 ≤ 𝜎22 If we consider 𝐸𝑗 to be the Mahalanobis circle M𝑗 (𝑆𝑗 , 𝑟𝑗 , 𝛴𝑗 ) with the
𝜗= 2
arctan 𝜎11 −𝜎22
∈⎨ (11)
⎪[𝜋∕2, 3𝜋∕4), 𝜎12 ≤ 0 & 𝜎11 < 𝜎22 center 𝑆𝑗 , the radius 𝑟𝑗 , and the covariance matrix 𝛴𝑗 , our problem is:
⎪[3𝜋∕4, 𝜋), 𝜎12 < 0 & 𝜎11 ≥ 𝜎22 .
⎩ ∑
𝑚
argmin 𝐹 (𝑆, 𝑟, 𝛴), 𝐹 (𝑆, 𝑟, 𝛴) = min D(M𝑗 (𝑆𝑗 , 𝑟𝑗 , 𝛴𝑗 ), 𝑎𝑖 ), (21)
𝑆∈𝛥𝑘 1≤𝑗≤𝑘
𝑖=1
Proof. The first statement of the lemma can be checked directly. In 𝑟∈[0,𝑅]𝑘 ,𝛴∈𝑀 𝑘
2
order to prove the reverse, let us rewrite M(𝑆, 𝑟; 𝛴) in the form
√ where 𝑀2 is the set of positive definite symmetric matrices of order
det 𝛴 two.
𝑟2
(𝑆− 𝑥)𝑇 𝛴 −1 (𝑆 − 𝑥) = 1, i.e. (12)
( 2 )−1 Notice that (21) is a GOP with 6𝑘 independent variables. As already
(𝑆 − 𝑥)𝑇 √ 𝑟 𝛴 (𝑆 − 𝑥) = 1. (13) mentioned, solving this kind of problem by using the DIRECT global
det 𝛴
( ) optimization algorithm would be very inefficient since the objective
𝑟2
By eigenvalue decomposition of the matrix √ 𝛴 the second part function is symmetric and has at least 𝑘! stationary points, and the
det 𝛴
of the lemma is proved. □ DIRECT algorithm should find all of them.
The most popular method for searching for a locally optimal par-
Similarly to the distance from a point to a circle, we can define
tition is the well-known 𝑘-means algorithm. It was proved that this
several different distances from a point 𝑎 to the ellipse, i.e., an M-circle
(see e.g. Gander et al., 1994; Grbić et al., 2016; Scitovski & Marošević, algorithm can find a globally optimal partition in the case when the
2014; Uteshev & Goncharova, 2018): initial approximation is very close to it. Bearing this in mind, we will
use the DIRECT algorithm only to find a good initial approximation
D1 (M(𝑆, 𝑟, 𝛴), 𝑎) = |‖𝑆 − 𝑎‖𝛴 − 𝑟| (LAD-distance), (14) and then use a modified 𝑘-means algorithm with M-circle-centers.
2
D2 (M(𝑆, 𝑟, 𝛴), 𝑎) = (‖𝑆 − 𝑎‖𝛴 − 𝑟) (TLS-distance), (15)
3.1.1. The adaptive mahalanobis 𝑘-closest M-circle-centers algorithm
D(M(𝑆, 𝑟, 𝛴), 𝑎) = (|‖𝑆 − 𝑎‖2𝛴 2 2
−𝑟 ) (algebraic distance). (16) A modification of the 𝑘-means algorithm for the case when the
The algebraic distance D is most frequently used in applications be- cluster centers are M-circles is called the Adaptive Mahalanobis 𝑘-
cause the minimizing function D is differentiable, and therefore it will closest M-circle-centers algorithm (KMCC); details can be found in (Grbić
also be used in our paper. et al., 2016; Marošević & Scitovski, 2015; Scitovski et al., 2021). The
Now we can write GOP (6) as follows: algorithm can be described in two steps which are repeated iteratively:
∑
𝑚
argmin D(M(𝑆, 𝑟, 𝛴), 𝑎𝑖 ), (17) Algorithm 1. (Adaptive Mahalanobis 𝑘-closest M-circle-centers
𝑆∈𝛥
𝑟∈[0,𝑅],𝛴∈𝑀2
𝑖=1 algorithm (KMCC))

3
R. Scitovski et al. Expert Systems With Applications 222 (2023) 119853

Step A: (Assignment step) Given a set  and 𝑘 different M-circles

M1 (𝑆1 , 𝑟1 , 𝛴1 ), … , M𝑘 (𝑆𝑘 , 𝑟𝑘 , 𝛴𝑘 ), apply the minimal distance
principle to define the clusters

𝜋𝑗 ∶= 𝜋𝑗 (M𝑗 ) = {𝑎 ∈  ∶ D(M𝑗 (𝑆𝑗 , 𝑟𝑗 , 𝛴𝑗 ), 𝑎)

≤ D(M𝑠 (𝑆𝑠 , 𝑟𝑠 , 𝛴𝑠 ), 𝑎), ∀𝑠 ≠ 𝑗}; (22)

Step B: (Update step) Given a partition 𝛱 = {𝜋1 , … , 𝜋𝑘 } of the set

, one can define the corresponding M-circle-centers M ̂𝑗
(𝑆̂𝑗 , 𝑟̂𝑗 , 𝛴̂ 𝑗 ), 𝑗 = 1, … , 𝑘 by solving the following GOP s:
∑
argmin D(M(𝑆𝑗 , 𝑟𝑗 , 𝛴𝑗 ), 𝑎), (23)
𝑆𝑗 ∈𝛥, 𝑟𝑗 ∈[0,𝑅], 𝛴𝑗 ∈𝑀2 𝑎∈𝜋𝑗

and calculate the objective function value 𝐹 given by (21). Fig. 1. Implementation of the KMCC algorithm.
̂ 𝑗 (𝑆̂𝑗 , 𝑟̂𝑗 , 𝛴̂ 𝑗 ), 𝑗 = 1, … , 𝑘.
Set M𝑗 (𝑆𝑗 , 𝑟𝑗 , 𝛴𝑗 ) ∶= M

Input: & Schnabel, 1996; Nocedal & Wright, 2006) since we are able to
 ⊂ 𝛥 (data set), 𝜖𝐵 > 0 determine a very favorable initial approximation similarly to what we
did in Section 2.3 by solving problem (17) for one ellipse.
Define M-circles Similarly to the classical 𝑘-means algorithm, it can be shown that
M0𝑗 (𝑆𝑗0 , 1, 𝐈), 𝑗 = 1, … , 𝑘 the sequence of the objective function values 𝐹𝑡 = 𝐹 (𝑆 (𝑡) , 𝑟(𝑡) , 𝛴 (𝑡) ),
𝑡 = 0, 1, … decreases monotonically (Späth, 1983). Hence, Algorithm
Set M𝑗 ∶= M0𝑗 , 𝑗 = 1, … , 𝑘, 1 can be stopped when the following condition is met for some small
and calculate 𝐹0 given by (21) 𝜖𝐵 > 0 (say 0.005) (see Bagirov et al., 2011):
𝐹𝑡−1 −𝐹𝑡
𝐹𝑡
< 𝜖𝐵 . (25)
For 𝑗 = 1, … , 𝑘 define 𝜋𝑗 ∶= 𝜋𝑗 (M𝑗 )
= {𝑎 ∈  ∶ D(M𝑗 (𝑆𝑗 , 𝑟𝑗 , 𝛴𝑗 ), 𝑎) Determining an initial approximation according to (24) and the
≤ D(M𝑠 (𝑆𝑠 , 𝑟𝑠 , 𝛴𝑠 ), 𝑎), ∀𝑠 ≠ 𝑗}
stopping criterion (25) was included in the given flowchart for the
KMCC algorithm.
̂ 𝑗 (𝑆̂𝑗 , 𝑟̂𝑗 , 𝛴̂ 𝑗 )
For each 𝜋𝑗 define M-circle-center M
Example 1. Let  be a data set coming from four ellipses shown
Set M𝑗 ∶= M ̂ 𝑗 , 𝑗 = 1, … , 𝑘, in Fig. 1a. This data set can be uploaded from DATAELL-1 (see Sec-
and calculate 𝐹1 given by (21) tion 4.1). By solving problem (24) we obtain four unit circles
(see Fig. 1a) used as an initial approximation for the KMCC algorithm.
After that, the KMCC algorithm gives the solution shown in Fig. 1b.
YES
(𝐹0 − 𝐹1 )∕𝐹1 > 𝜖𝐵 Set 𝐹0 = 𝐹1

NO 3.2. The number of ellipses is not known in advance

Output:
𝐹1 , 𝛱 ∶= {𝜋1 , … , 𝜋𝑘 }, {M1 , … , M𝑘 } For the given data point set  by using the Incremental method
(see e.g. Bagirov et al., 2011; Scitovski & Sabo, 2019) we are able
Algorithm 1 (Adaptive Mahalanobis k-closest M-circle-centers algorithm to successfully determine optimal partitions with 1, 2, … clusters for
(KMCC)) the case when cluster-centers are ellipses (M-circles). Among them, we
need to determine the partition with the most appropriate number of
clusters, the most appropriate partition (MAPart). In this paper, we will
The algorithm starts with Step A, for which a good initial approxi- define a new index specialized in finding a MAPart for the case when
mation for M-circle-centers will be necessary. For that purpose, we will cluster-centers are ellipses (M-circles).
search for an optimal partition with unit circles as cluster-centers, i.e., Let 𝛱 ⋆ = {𝜋1⋆ , … , 𝜋𝑘⋆ } be the optimal partition of the set  =
we will solve the following simple GOP {𝑎 = (𝑥𝑖 , 𝑦𝑖 )𝑇 ∶ 𝑖 = 1, … , 𝑚} ⊂ R2 with ellipses 𝐸𝑗⋆ = 𝐸𝑗⋆ (𝑆𝑗 , 𝜉𝑗 , 𝜂𝑗 , 𝜗𝑗 ),
𝑖

∑
𝑚 𝑗 = 1, … , 𝑘 as cluster-centers.
𝑆 0 ∈ argmin 𝛷(𝑆), 𝛷(𝑆) = min D(𝑗 (𝑆𝑗 , 1), 𝑎𝑖 ), (24) For each cluster 𝜋𝑗⋆ we know its number of elements |𝜋𝑗⋆ | and the
𝑆∈𝛥𝑘 1≤𝑗≤𝑘
𝑖=1
corresponding center, i.e., its M-circle-center 𝐸𝑗⋆ . Since 𝐸𝑗⋆ is an ellipse
with D(𝑗 (𝑆𝑗 , 1), 𝑎𝑖 ) = (‖𝑆𝑗 − 𝑎𝑖 ‖2 − 1)2 being the algebraic distance with semi-axes 𝜉 and 𝜂, its length (circumference) is
between an ordinary unit circle 𝑗 with the center 𝑆𝑗 and the point √
𝜋∕2
𝑎𝑖 ∈ . GOP (24) is a simpler problem that has 2𝑘 variables and can be |𝐸𝑗⋆ | = 4 𝜉𝑗2 cos2 𝑡 + 𝜂𝑗2 sin2 𝑡 𝑑𝑡. (26)
solved by using the DIRECT algorithm. Since the solution of GOP (24) ∫0
is used only as an initial approximation in the KMCC algorithm, it is This is an elliptic integral of the second kind and cannot be solved ele-
enough to run only a few iterations of the DIRECT algorithm. mentarily, but a very good approximation is the well-known Ramanujan
After that, we can form an M-circle M0𝑗 (𝑆𝑗0 , 1, 𝐈), 𝑗 = 1, … , 𝑘 as the approximation (see also Sabo & Scitovski, 2020):
initial approximation in Step A, where 𝐈 is an identity matrix of order
( 3ℎ ) (𝜉𝑗 −𝜂𝑗 )2
two. Notice that in this way the centers of searched ellipses will be well |𝐸𝑗⋆ | ≈ 𝜋(𝜉𝑗 + 𝜂𝑗 ) 1 + √ , ℎ= (𝜉𝑗 +𝜂𝑗 )2
. (27)
10+ 4−3ℎ
positioned.
Solutions to problem (23) can be found by using some local opti- Let 𝛱 ⋆ = {𝜋1⋆ , … , 𝜋𝑘⋆ } be a globally optimal partition with M-circle-
mization method as Nelder–Meade or Quasi-Newton (see e.g. Dennis centers 𝐸1⋆ , … , 𝐸𝑘⋆ . For each cluster 𝜋𝑗⋆ we define:

4
R. Scitovski et al. Expert Systems With Applications 222 (2023) 119853

Fig. 2. Four selected examples from DATAELL-1 with 𝑘 = 2, 3, 4, 5 ellipses.

Fig. 3. Four selected examples from DATAELL-2 with 𝑘 = 2, 3, 4, 5 ellipses.

• the mean squared deviation of its points from M-circle-center 𝐸𝑗⋆ Each of 100 sets in DATAELL-1 is constructed as follows. First,
we randomly choose an integer 𝑘 ∈ {2, 3, 4, 5}. After that, we choose
1
∑ 𝑘 random ellipse-centers 𝑆1 , … , 𝑆𝑘 ∈ [1.5, 8.5]2 ⊂ [0, 10]2 that are
𝖣𝖾𝗏𝗂𝖺𝗍𝗂𝗈𝗇𝗌(𝜋𝑗⋆ ) ∶= |𝜋𝑗⋆ |
D(𝐸𝑗⋆ , 𝑎), (28)
mutually distant by at least 3. Semi-axes 𝜉𝑗 , 𝜂𝑗 > 0 of ellipses will be
𝑎∈𝜋𝑗⋆
taken from  (0.5, 2), and the angle of rotation 𝜗𝑗 ∈  (0, 𝜋). Ellipses
where D(𝐸𝑗⋆ , 𝑎) is algebraic distance (16) from the M-circle-center 𝐸𝑗 (𝐶𝑗 , 𝜉𝑗 , 𝜂𝑗 , 𝜗𝑗 ), 𝑗 = 1, … , 𝑘 constructed in this way will probably
𝐸𝑗⋆ to 𝑎 ∈ 𝜋𝑗⋆ ; intersect a bit.
• the density around the M-circle-center 𝐸𝑗⋆ Then, we randomly choose integers 𝑚𝑗 ∈ {180, … , 220} uniformly
distributed points (see Grbić et al., 2016) on the ellipse 𝐸𝑗 and at each
|𝜋𝑗⋆ |
𝖣𝖾𝗇𝗌𝗂𝗍𝗒(𝜋𝑗⋆ ) ∶= |𝐸𝑗⋆ |
. (29) point in the direction of the normal vector we take a random point
𝑎𝑖 ∈  (0, 𝜎 2 ), 𝜎 2 = 0.05. Fig. 2 shows four such sets.
Due to the mutual intersection of ellipses, these properties can differ DATAELL-1 is loaded in Mathematica as:
slightly.
The closer a partition 𝛱 ⋆ to the MAPart the smaller 𝖣𝖾𝗏𝗂𝖺𝗍𝗂𝗈𝗇𝗌(𝜋𝑗⋆ )
In[1]:= SetDirectory[NotebookDirectory[]];
of all 𝜋𝑗⋆ and the greater 𝖣𝖾𝗇𝗌𝗂𝗍𝗒(𝜋𝑗⋆ ). Therefore it makes sense to define
DATAELL = ReadList["DATAELL-1.txt"];
a new index specialized in finding an optimal MAPart with M-circles
as center-centers:
For the selected 𝑖𝑖∈ {1, … , 100}, DATAELL[[ii,1]] and
𝑘 𝖣𝖾𝗏𝗂𝖺𝗍𝗂𝗈𝗇𝗌(𝜋 ⋆ )
∑ ∑
𝑘 ∑
𝑗 |𝐸𝑗⋆ | DATAELL[[ii,2]] represent the data point set and the correspond-
⋆
𝙶𝙾(𝑘) ∶= 𝙶𝙾(𝛱 (𝑘)) = 𝑘 =𝑘 2 D(𝐸𝑗⋆ , 𝑎).
𝖣𝖾𝗇𝗌𝗂𝗍𝗒(𝜋𝑗⋆ ) |𝜋𝑗⋆ | ing ellipse parameters (centers, semi-axes, angles of rotation), respec-
𝑗=1 𝑗=1 𝑎∈𝜋𝑗⋆
tively.
(30) Data set DATAELL-2 is constructed in a similar way. 𝑘 random
ellipse-centers will be chosen in the square [1.5, 8.5]2 ⊂ [0, 10]2 , within
This index will be called the Geometrical Objects-index (GO-
a range of at least 2.5. Semi-axes 𝜉, 𝜂 > 0, will be taken from  (.5, 2.5),
index). It is clear that a small value of the GO-index indicates the
while the angle of rotation 𝜗 remains from  (0, 𝜋). A few of such
optimal MAPart. Its main advantage over other known indexes is that
ellipses are shown in Fig. 3.
it is especially used for searching for a MAPart with clusters whose
Ellipse detection from data set DATAELL-1 or DATAELL-2 can be
centers can be various geometrical objects (circles, ellipses, lines, etc.).
done by our Mathematica software available on: http://clusters.mathos.
4. Numerical experiments unios.hr/
The proposed method has also been tested on several data sets
4.1. Data sets and our publicly available software coming from real-world images. The corresponding algorithm was im-
plemented in C++ by using two open-source libraries. In the case of
The proposed method will be tested on 100 artificial data sets com- recognizing ellipses in real-world images, first we need to preprocess
ing from preponderantly well separated ellipses (DATAELL-1), which the edges as the input for the algorithms are data points. Edge detection
can be downloaded from: http://clusters.mathos.unios.hr/modules/ is done by using the Canny filter from OpenCV (Bradski, 2000). For
DATAELL-1.txt, and 100 more complex artificial data sets coming solving the corresponding optimization problems we used the NLopt
from mutually intersecting ellipses (DATAELL-2), which can be down- library (Johnson, 2011) for nonlinear optimization, providing an in-
loaded from: http://clusters.mathos.unios.hr/modules/DATAELL-2.txt, terface for different optimization routines. We use its implementation
and from data sets coming from real-world images. for the DIRECT algorithm (Jones et al., 1993) and the Nelder–Mead

5
R. Scitovski et al. Expert Systems With Applications 222 (2023) 119853

algorithm (Nelder & Mead, 1965). A demo version of our software is The algorithm could speed up significantly if we had an index
available on: http://cs.mathos.unios.hr/~pnikic/ells/. that decreases monotonically to the optimal partition. The Davies–
Bouldin index (see e.g. Vendramin et al., 2009) does not have this
4.2. Searching for a MAPart property, and a new GO-index almost always exemplifies this property.
If the index with a monotonically decreasing property were used, new
As the objective function value (6) decreases by increasing the partitions with one cluster more would be calculated as long as the
number of clusters, starting with a 1-partition, by applying the KMCC al- index value decreased. This would speed up the algorithm. Since it
gorithm we calculate new partitions with one cluster more until the cannot be guaranteed that the GO-index will always show the property
relative objective function value becomes smaller than 𝜖 > 0 given in of monotonicity, in that case we should take a risk that it may happen
advance. This is a well-known incremental algorithm (see e.g. Bagirov that we find only a locally optimal partition. Algorithm 3 is constructed
et al., 2011; Scitovski & Sabo, 2019). with the assumption of monotonicity of the GO-index.
For each partition obtained, we can determine the values of one Notice that the KMCC algorithm used in Algorithm 2 and Algo-
of the unknown indexes, e.g., Davies–Bouldin, Calinski–Harabasz (Liu rithm 3 begins with determining an initial approximation obtained by
et al., 2013; Vendramin et al., 2009) and try to determine the solving GOP (24).
MAPartbased on them. In numerous examples, it is shown that for The corresponding flowchart is given along with Algorithm 3.
this purpose it is better to use the proposed GO-index as it is used
especially for searching for the MAPart with clusters whose centers
Algorithm 3 (Searching for an optimal MAPart)
can be various geometrical objects.
Algorithm 2 is an incremental algorithm, in each step of which we Input:  ⊂ 𝛥 {Data set};
calculate the value of the GO-index. The smallest value of the GO-index 1: Set 𝑘 = 1 and call the KMCC algorithm;
points to the MAPart. 2: Set 𝑘 = 𝑘 + 1 and call the KMCC algorithm;
The corresponding flowchart for Algorithm 2 is given below. 3: while (GO(𝑘) < GO(𝑘 − 1) do
4: Set 𝑘 = 𝑘 + 1 and call the KMCC algorithm;
5: end while
Algorithm 2 (Searching for an optimal MAPart) Output: {𝑘, 𝛱 ⋆ , {𝐸1⋆ , … , 𝐸𝑘⋆ }}.
Input:  ⊂ 𝛥 {Data set}; 𝜖 > 0;
1: Set 𝑡 = 1 and call the KMCC algorithm;
Input:
2: Set 𝑡 = 𝑡 + 1 and call the KMCC algorithm;  ⊂ 𝛥 (data set)
3: if (𝐹𝑡−1 − 𝐹𝑡 )∕𝐹1 > 𝜖, then
4: Set 𝑡 = 𝑡 + 1 and call the KMCC algorithm and go to Step 3; Set 𝑘 = 1
5: end if and call KMCC algorithm
6: Find a 𝑘-partition 𝛱 ⋆ (𝑘 ≤ 𝑡) with the smallest GO-index and denote
the corresponding cluster-centers by 𝐸1⋆ , … , 𝐸𝑘⋆ ; Set 𝑘 = 𝑘 + 1
Output: {𝑘, 𝛱 ⋆ , {𝐸1⋆ , … , 𝐸𝑘⋆ }}. and call KMCC algorithm

YES
Input: GO(𝑘) < GO(𝑘 − 1)
 ⊂ 𝛥 (data set), 𝜖 > 0
NO
Set 𝑡 = 1 Output:
and call KMCC algorithm 𝑘, 𝛱 ⋆ , {𝐸1⋆ , … , 𝐸𝑘⋆ }

Set 𝑡 = 𝑡 + 1
Algorithm 3 (Searching for an optimal MAPart)
and call KMCC algorithm
4.3. Testing the method and the algorithms
YES
(𝐹𝑡−1 − 𝐹𝑡 )∕𝐹1 > 𝜖 The method for solving the multiple ellipse detection problem given
in Algorithm 3 will be tested1 first on 100 sets of data points from
NO DATAELL-1 (see Section 4.1). These are simple examples of ellipses
Find a 𝑘-partition, 𝑘 < 𝑡, with the smallest GO-index
which do not overlap or intersect just a little bit. Table 1 shows the
and denote corresponding cluster-centers 𝐸1⋆ , … , 𝐸𝑘⋆ results.
Let us recall that Algorithm 3 is constructed with the assumption
of monotonicity of the GO-index. Table 1 shows the number of ellipses
Output:
𝑘, 𝛱 ⋆ , {𝐸1⋆ , … , 𝐸𝑘⋆ }
detected by Algorithm 3, especially for data sets coming from 2, 3, 4,
or 5 ellipses. As can be seen, Algorithm 3 recognizes very well ellipses
the set of data points  was generated from. Necessary CPU-time is
Algorithm 2 (Searching for an optimal MAPart) very low and consists of the CPU-time necessary to run 15 iterations of
the DIRECT algorithm (for obtaining an initial approximation) and the
Example 2. For the data set from Example 1 coming from four ellipses CPU-time necessary to run the KMCC algorithm. Monotonicity of the
shown in Fig. 2c, by applying Algorithm 2 we determine optimal 𝑘- GO-index contributed significantly to this phenomenon. On the other
partitions for 𝑘 = 1, 2, 3, 4, 5, and calculate the objective function values hand, this was the reason why in a small number of cases the algorithm
𝐹𝑘 , 𝑘 = 2, 3, 4, 5 (according to (6)), and the value 𝙶𝙾(𝑘) of the GO-index. has not recognized the ellipses searched for.
The smallest value of the GO-index is obtained for 𝑘 = 4. This means
that the GO-index points to the 4-partition as the most appropriate one.
Comparing the obtained M-circle-centers with the original one, we can 1
All evaluations were performed on a computer with a 2.50 GHz
see that we have found a globally optimal partition (see Fig. 4). Intel(R) Core(TM) i5-2520M CPU with 4GB of RAM.

6
R. Scitovski et al. Expert Systems With Applications 222 (2023) 119853

Fig. 4. Implementation of Algorithm 2.

Table 1
Efficiency of Algorithm 3 in the case of separated ellipses.
No. of ellipses Ellipses detected CPU-time
0 1 2 3 4 5 6 (DIRECT) (k-means) Total
2 0 0 20 0 0 0 0 0.01166 0.03388 0.04554
3 0 0 0 25 0 0 0 0.03963 0.08338 0.12301
4 1 0 0 0 26 0 0 0.09531 0.11201 0.20732
5 5 1 0 1 0 21 0 0.09542 0.15709 0.25251

Table 2
Efficiency of Algorithm 2 in the case of separated ellipses.
No. of ellipses Ellipses detected CPU-time
0 1 2 3 4 5 6 (DIRECT) (k-means) Total
2 0 0 20 0 0 0 0 0.01544 0.03781 0.05325
3 0 0 0 25 0 0 0 0.03723 0.08193 0.11916
4 0 0 0 0 27 0 0 0.09917 0.11905 0.21822
5 0 0 0 1 0 27 0 0.11779 0.17902 0.29681

Table 3
Efficiency of Algorithm 3 in the case of intersecting ellipses.
No. of ellipses Ellipses detected CPU-time
0 1 2 3 4 5 6 (DIRECT) (k-means) Total
2 0 1 19 0 0 0 0 0.01221 0.06475 0.07696
3 0 1 3 21 0 0 0 0.04443 0.15606 0.20049
4 8 2 2 2 13 0 0 0.07108 0.19019 0.26127
5 5 3 3 8 1 8 0 0.09005 0.25748 0.34753

If we used Algorithm 2 instead of Algorithm 3 (i.e., if we abandon set  and the algorithm should be run again with the rest of the data
the application of the monotonicity property of the GO-index) in these on a new, simpler set.
examples, recognition results would be even better, i.e., there would
Finally, let us say a few words about a comparison of our method
be only one example with 5 ellipses in which our algorithm recognizes
with some other known algorithms. As shown in Scitovski and Maroše-
3 out of 5 ellipses, but necessary CPU-time slightly increases (see Ta-
vić (2014), for solving such problems, the method based on the Hough
ble 2). It can be said that in 99% of the cases our algorithm recognizes
transformation is inferior in relation to the incremental algorithm,
such ellipses. This result also shows high efficiency of the GO-index.
and the incremental algorithm performed worse in relation to the
Algorithm 3 will also be tested on more complex examples of
statistically based method (Grbić et al., 2016). In similar examples, the
mutually intersecting ellipses. Table 3 shows the results. As expected,
new proposed method shows significant improvements.
in this case, the number of successfully detected ellipses is smaller, and
The well-known EDCircles algorithm, proposed by Akinlar and
necessary CPU-time has not increased significantly (see Table 3).
Topal (2013), has very good characteristics for the case of recognizing
This example points to the possibility of reconstructing Algorithm 2,
ellipses with clear edges, but it does not recognize ellipses with unclear
i.e., Algorithm 3, in order to increase the number of detected ellipses.
or noisy edges. This is the advantage of the new proposed method.

Remark 1. The efficiency of Algorithm 2 and Algorithm 3 could be

4.4. Recognizing ellipses in real-world images
improved. When searching for an optimal 𝑘-partition and determining
the corresponding GO-index for every cluster 𝜋𝑗⋆ , by means of (29)–(27)
The use of the new proposed method will also be illustrated by a
we can determine the following indicator:
few examples of ellipse recognition in real-world images. Other images
|𝐸𝑗⋆ | ∑ can be taken from the Caltech 256 data set (Griffin et al., 2007). For
2 D(𝐸𝑗⋆ , 𝑎). (31)
|𝜋𝑗⋆ | that purpose, it is necessary to perform edge detection preprocessing
𝑎∈𝜋𝑗⋆
by using the Canny filter from OpenCV (Bradski, 2000). The original
A small value of number (31) will indicate if the ellipse-center 𝐸𝑗⋆ images and detected ellipses are shown in Figs. 5 and 6.
has been recognized as one of the ellipses searched for. In that case, Table 4 gives the obtained results: the number of detected ellipses
similarly to Grbić et al. (2016), cluster 𝜋𝑗⋆ can be taken out from the and CPU-time for two examples of real-world images shown in Figs. 5

7
R. Scitovski et al. Expert Systems With Applications 222 (2023) 119853

Fig. 5. Real-world image detection.

Fig. 6. Real-world image detection.

Table 4 good initial approximation. For that purpose, a solution to the auxiliary
Efficiency of Algorithm 2 in real-world images.
problem is searched for, i.e., the multiple unit circles detection problem
Ellipses CPU-time (see Scitovski & Sabo, 2019), and for obtaining a good initial approxi-
detected
mation it suffices to run only a few iterations of the global optimization
(DIRECT) (k-means) Total
algorithm DIRECT.
Coins (Fig. 5) 6 0.01544 0.03781 0.05325
Cups (Fig. 6) 6 0.03723 0.08193 0.11916 In the case when the number of ellipses is not known in advance, for
the data points  the MAPart is searched for by using the incremental
algorithm (see e.g. Bagirov et al., 2011; Scitovski et al., 2021) and the
new proposed Geometrical Objects-index (GO-index).
and 6. Calculations were done by using Algorithm 2 implemented in
C++. Numerical experiments show high efficiency of the proposed
method, and the new GO-index that identifies the MAPart proved to
be significantly more efficient than other indexes. Unlike numerous
5. Conclusions and future work
algorithms presented in the literature, which recognize ellipses with
clear edges well, but do not recognize ellipses with unclear or noisy
A new method for solving the multiple ellipse detection problem edges, our algorithm is specifically used for such situations.
is proposed in the paper based on the data point set  coming from
a number of ellipses in the plane not known in advance. In so doing, The method proposed in this paper can be applied to multiple
data points are usually contaminated with some noisy errors. ellipses detection in real-world images. It is thereby necessary to do
image preprocessing first, which is based on some edge detection
The method is based on center-based clustering whereby cluster-
method (see e.g. Mathematica computation system (Wolfram, 2020) or
centers are Mahalanobis circles. In this way, searching for an unknown
the Canny filter from OpenCV Bradski, 2000). Furthermore, our method
ellipses is reduced to the problem of searching for the most appropriate
proved to be very efficient in numerous real-world images, except in
partition with Mahalanobis circle-centers.
the case when the preprocessed image contains, in addition to ellipses,
In the case when the number of ellipses is known in advance, a
a large number of non-ellipsoidal geometrical objects.
solution to the multiple ellipse detection problem is searched for by
modifying the 𝑘-means algorithm for the case of Mahalanobis circles It is this shortcoming of the proposed method that we will try
as cluster centers. In this case, it is extremely important to have a to eliminate in our future work. Additional analysis and a possible

8
R. Scitovski et al. Expert Systems With Applications 222 (2023) 119853

correction of the proposed GO-index could also be part of our future Jones, D. R., & Martins, J. R. R. A. (2020). The DIRECT Algorithm–25 years later.
work. Journal of Global Optimization, 79, 521–566.
Jones, D. R., Perttunen, C. D., & Stuckman, B. E. (1993). Lipschitzian optimization
without the Lipschitz constant. Journal of Optimization Theory and Applications, 79,
CRediT authorship contribution statement 157–181.
Kashiha, M. A., Bahr, C., Ott, S., Moons, C. P., Niewold, T. A., & Berckmans, F. T.
Rudolf Scitovski: Conceptualization, Methodology, Software, Vali- D. (2014). Automatic monitoring of pig locomotion using image analysis. Livestock
dation, Formal analysis, Investigation, Data curation, Writing – review Science, 159, 141–148.
Kwon, B. K., Teng, Z., Roh, T. J., & Kang, D. J. (2016). Fast ellipse detection based on
& editing, Visualization, Supervision, Project administration. Kristian three point algorithm with edge angle information. International Journal of Control,
Sabo: Conceptualization, Methodology, Software, Validation, Formal Automation and Systems, 14, 804–813.
analysis, Investigation, Resources, Writing – review & editing, Vi- Liu, C., Chen, R., Chen, K., & Xu, J. (2022). Ellipse detection using the edges extracted
sualization, Supervision, Project administration, Funding acquisition. by deep learning. Machine Vision and Applications, 33(63).
Liu, Y., Li, Z., Xiong, H., Gao, X., Wu, J., & Wu, S. (2013). Understanding and enhance-
Patrick Nikić: Methodology, Software, Formal analysis, Investigation,
ment of internal clustering validation measures. IEEE Transaction on Cybernetics, 43,
Data curation, Visualization. Snježana Majstorović: Methodology, 982–994.
Validation, Writing – review & editing. Liu, Y. S., & Ramani, K. (2009). Robust principal axes determination for point-based
shapes using least median of squares. Computer-Aided Design, 41, 293–305.
Declaration of competing interest Liu, Y., Xie, Z., & Liu, H. (2019). Fast and robust ellipse detector based on edge
following method. IET Image Processing, 13, 2409–2419.
Lu, C., Xia, S., Shao, M., & Fu, Y. (2020). Arc-support line segments revisited: An
The authors declare that they have no known competing finan- efficient high-quality ellipse detection. IEEE Transactions on Image Processing, 29,
cial interests or personal relationships that could have appeared to 768–781.
influence the work reported in this paper. Maalek, R., & Lichti, D. D. (2021). Robust detection of non-overlapping ellipses
from points with applications to circular target extraction in images and cylinder
detection in point clouds. ISPRS Journal of Photogrammetry and Remote Sensing, 178,
Data availability 83–108.
Marošević, T., & Scitovski, R. (2015). Multiple ellipse fitting by center-based clustering.
Data will be made available on request. Croatian Operational Research Review, 6, 43–53.
Moshtaghi, M., Havens, T. C., Bezdek, J. C., Park, L., Leckie, C., Rajasegarar, S.,
Keller, J. M., & Palaniswami, M. (2011). Clustering ellipses for anomaly detection.
Acknowledgments
Pattern Recognition, 44, 55–69.
Nelder, J. A., & Mead, R. (1965). A simplex method for function minimization. The
The authors would like to thank the referees and the journal editors Computer Journal, 7, 308–313.
for their careful reading of the paper and insightful comments that Nocedal, J., & Wright, S. J. (2006). Numerical optimization (2nd ed.). Springer.
Prasad, D. K., Leung, M. K. H., & Quek, C. (2013). Ellifit: An unconstrained, non-
helped us improve the paper. Especially, the authors would like to
iterative, least squares based geometric ellipse fitting method. Pattern Recognition,
thank Mrs. Ivanka Ferčec for significantly improving the use of English 46, 1449–1465.
in the paper, and Prof. Dr Šime Ungar for drawing flowcharts of our Rueda, S., Fathima, S., Knight, C. L., Yaqub, M., Papageorghiou, A. T., Rahmatullah, B.,
algorithms by using the LATEX package tikz. Foi, A., Maggioni, M., Pepe, A., Tohka, J., Stebbing, R. V., McManigle, J. E.,
Ciurte, A., Bresson, X., Cuadra, M. B., Sun, C., Ponomarev, G. V., Gelfand, M.
S., Kazanov, M. D., .... Noble, J. A. (2013). Evaluation and comparison of current
References
fetal ultrasound image segmentation methods for biometric measurements: A grand
challenge. IEEE Transactions on Medical Imaging, 10, 1–16.
Akinlar, C., & Topal, C. (2013). Edcircles: A real-time circle detector with a false Sabo, K., Grahovac, D., & Scitovski, R. (2020). Incremental method for multiple line
detection control. Pattern Recognition, 46, 725–740. detection problem - Iterative reweighted approach. Mathematics and Computers in
Bagirov, A. M., Ugon, J., & Webb, D. (2011). Fast modified global 𝑘-means algorithm Simulation, 178, 588–602.
for incremental cluster construction. Pattern Recognition, 44, 866–876. Sabo, K., & Scitovski, R. (2020). Multiple ellipse detection by using RANSAC and
Bradski, G. (2000). The openCV library. Dr. Dobb’s Journal of Software Tools. DBSCAN method. In Proceedings of the 9th international conference on pattern
Dennis, J. J., & Schnabel, R. (1996). Numerical methods for unconstrained optimization recognition applications and methods (pp. 129–135). http://dx.doi.org/10.5220/
and nonlinear equations. Philadelphia: SIAM. 0008879301290135.
Dong, H., Prasad, D. K., & Chen, I. (2018). Accurate detection of ellipses with false Scitovski, R., & Marošević, T. (2014). Multiple circle detection based on center-based
detection control at video rates using a gradient analysis. Pattern Recognition, 81, clustering. Pattern Recognition Letters, 52, 9–16.
112–130. Scitovski, R., & Sabo, K. (2019). Application of the DIRECT algorithm to searching
Fitzgibbon, A., Pilu, M., & Fisher, R. B. (1999). Direct least square fitting of ellipses. for an optimal 𝑘-partition of the set A and its application to the multiple circle
IEEE Transactions on Pattern Analysis and Machine Intelligence, 21, 476–480. detection problem. Journal of Global Optimization, 74, 63–77.
Fornaciari, M., Prati, A., & Cucchiara, R. (2014). A fast and effective ellipse detector Scitovski, R., Sabo, K., Martínez-Álvarez, F., & Ungar, Š. (2021). Cluster analysis and
for embedded vision applications. Pattern Recognition, 47, 3693–3708. applications. Springer.
Gander, W., Golub, G. H., & Strebel, R. (1994). Least-squares fitting of circles and Späth, H. (1983). Cluster-formation und analyse. München: R. Oldenburg Verlag.
ellipses. BIT, 34, 558–578. Uteshev, A. Y., & Goncharova, M. V. (2018). Point-to-ellipse and point-to-ellipsoid
Grbić, R., Grahovac, D., & Scitovski, R. (2016). A method for solving the multiple distance equation analysis. Journal of Computational and Applied Mathematics, 328,
ellipses detection problem. Pattern Recognition, 60, 824–834. 232–251.
Grbić, R., Nyarko, E. K., & Scitovski, R. (2013). A modification of the DIRECT Vendramin, L., Campello, R. J. G. B., & Hruschka, E. R. (2009). On the comparison of
method for Lipschitz global optimization for a symmetric function. Journal of Global relative clustering validity criteria. In Proceedings of the SIAM international conference
Optimization, 57, 1193–1212. on data mining (pp. 733–744). SIAM.
Griffin, G., Holub, A., & Perona, P. (2007). Caltech-256 object category database: Technical Wang, G., Ren, G., Wu, Z., Zhao, Y., & Jiang, L. (2014). A fast and robust ellipse-
report, Caltech, Available http://authors.library.caltech.edu/7694S. detection method based on sorted merging. The Scientific World Journal, 2014, ID
Jia, Q., Fan, X., Luom, Z., Song, L., & Qiu, T. (2017). A fast ellipse detector using 481312.
projective invariant pruning. IEEE Transactions on Image Processing, 26, 3665–3679. Wolfram, S. (2020). An elementary introduction to the Wolfram language (2nd ed.).
Johnson, S. G. (2011). The NLopt nonlinear-optimization package. http://ab-initio.mit. Champaign, Illinois: Wolfram Research, Inc., Version 12.0.
edu/nlopt. Wynn, T. J., & Stewart, S. A. (2005). Comparative testing of ellipse-fitting algorithms:
Implications for analysis of strain and curvature. Journal of Structural Geology, 27,
1973–1985.