#  Rosetta Stone at the Arabic Reverse Dictionary Shared Task:
A Hop From Language Modeling To Word–Definition Alignment

Ahmed ElBakry*
Microsoft
<EMAIL_ADDRESS>
&Mohamed Gabr
Microsoft
<EMAIL_ADDRESS>
Muhammad ElNokrashy
Microsoft
<EMAIL_ADDRESS>
&Badr AlKhamissi*
EPFL
<EMAIL_ADDRESS>
###### Abstract

A Reverse Dictionary is a tool enabling users to discover a word based on its
provided definition, meaning, or description. Such a technique proves valuable
in various scenarios, aiding language learners who possess a description of a
word without its identity, and benefiting writers seeking precise terminology.
These scenarios often encapsulate what is referred to as the "Tip-of-the-
Tongue" (TOT) phenomena. In this work, we present our winning solution for the
Arabic Reverse Dictionary shared task. This task focuses on deriving a vector
representation of an Arabic word from its accompanying description. The shared
task encompasses two distinct subtasks: the first involves an Arabic
definition as input, while the second employs an English definition. For the
first subtask, our approach relies on an ensemble of finetuned Arabic BERT-
based models, predicting the word embedding for a given definition. The final
representation is obtained through averaging the output embeddings from each
model within the ensemble. In contrast, the most effective solution for the
second subtask involves translating the English test definitions into Arabic
and applying them to the finetuned models originally trained for the first
subtask. This straightforward method achieves the highest score across both
subtasks.111https://github.com/bkhmsi/RashidRevDict

**footnotetext: Equal contribution

##  1 Introduction

The Tip-of-the-Tongue phenomena, as explained by the authors of Brown and
McNeill (1966), is “a state in which one cannot quite recall a familiar word
but can recall words of similar form and meaning”. A straightforward way to
solve this problem, is to have a reverse dictionary; a system that takes a
description provided by the user as an input, and outputs the word (Bilac et
al., 2004).

The initial solutions were heuristic-based. In their work, Shaw et al. (2013)
suggested a method where the tokens in the user-provided description are
compared to all dictionary definitions. The system then returns the word with
the highest token match. Their method implements different retrieval
efficiency tweaks to overcome the issue of excessive time complexity resulting
from the comparison operation.

Recent approaches employ neural-based models, since that are capable of better
capturing the semantics of an input description, in contrast to the earlier
solutions mentioned, which relied on word overlap. In their work, Hill et al.
(2015) suggest utilizing a recurrent-neural-network (RNN) to generate a vector
representation based on the provided definition. This representation is then
compared against a set of word embeddings to select the closest word to return
to the user.

The issue of low-frequency words is one of the main challenges of building a
reverse dictionary, since these words are the ones that are the less trained
and thus have a worse representation compared to more frequent words. Zhang et
al. (2019) tackle this problem by handcrafting predictors that extract
features inspired by the thought process undergone by humans to get a word
given its description.

Polysemy, which is the coexistence of many possible meanings for a word, is
another obstacle when building reverse dictionaries. Most of the previously
mentioned solutions rely on a a set of static word representation builders
such as Word2Vec, which hinders the accuracy of such models. This motivates
the use of pretrained language models to produce embeddings that vary based on
context. The authors of Yan et al. (2020) probed BERT Devlin et al. (2019) to
predict the word representation, alleviating the issue of polysemy.

| Example Word | Example Definition
---|---|---
Task 1 |  t.hnnana ‘lyh |  tar.ham, ta‘.taf ‘lyh w ra.hmah
Task 2 |  zwwar alkalaama | To knowingly and willfully make a false statement of witness while in court
Table 1: Word-Definition Pairs Illustrating Subtasks 1 and 2.

Reverse dictionaries can also be cross-lingual; where one aims to retrieve a
word in language X based on a description provided in language Y. Employing
any of the previously mentioned solutions for a multilingual context
necessitates the alignment of word vectors across different languages, a
challenging task even for two languages, not to mention when dealing with
multiple languages. The authors of Chen et al. (2018) built a collection of
bilingual reverse dictionaries using Wiktionary. Other solutions used existing
multilingual models, such as mBERT, to reduce the issue of cross-lingual
alignment Yan et al. (2020).

The shared task of the Arabic Reverse Dictionary provides a set of words,
along with their  AySGNS  Mikolov et al. (2013) and  AyELECTRA  Clark et al.
(2020) vector representations, and their corresponding definition, in both
Arabic and English. A set of Arabic-English word mappings is also supplied to
help in building an alignment scheme. The goal of subtask 1 is to predict the
AySGNS  and  AyELECTRA  embeddings of the set of Arabic words, given the input
Arabic definition. Subtask 2 has the same goal except that an English
definition is provided instead of Arabic.

The shared task setup poses multiple obstacles that our solutions attempt to
overcome: (1) the small size of the set of aligned words, (2) the blackbox
nature of the  AySGNS  and  AyELECTRA  word embedding generation pipeline.

Our solution simply finetunes multiple Arabic BERT-based pretrained models to
predict an embedding for each word.

##  2 Datasets

The provided data can be categorised into three distinct datasets.

  1. 1.

The Arabic Language Dictionary is a dataset with $58,010$ entries, where each
of datapoint contains a word, an  AyELECTRA  embedding, an  AySGNS  embedding,
a  Aygloss  (definition of the word), a  AyPOS  tag, an  AyID  and an
AyEnglish ID  where applicable to link with the alignment data.

  2. 2.

The English Language Dictionary dataset has $63,596$ datapoints, with the same
columns as the Arabic Dictionary except that the embeddings are obtained from
English words and not Arabic.

  3. 3.

The English Arabic Mapped Dictionary has $4,355$ datapoints in total. Each
point has the Arabic and English glosses, Arabic and English IDs, Arabic and
English words, and the Arabic embeddings.

The first and third datasets are split into  Aytrain ,  Aydev  and  Aytest
sets by the organizers. The English language dictionary however, isn’t
provided with such divisions. Therefore, we manually split the English
dictionary ourselves. Table 2 shows the split sizes of each dataset. The
English dictionary was divided into two sets only,  Aytrain  and  Aydev ,
since there was no need for a test set in our case, and no submission to be
made with this dictionary.

| Train | Dev | Test
---|---|---|---
Ar Dict | 45,200 | 6,400 | 6,410
Ar-En Map | 2,843 | 299 | 1,213
Ar Dict | 50,877 | 12,719 | N/A
Table 2: Statistics about Data Sizes

##  3 System

###  3.1 Subtask 1: Arabic Definitions to Arabic Embeddings

In this subtask, we finetune four Arabic BERT-based pretrained models. Namely:
(1) MARBERTv2 (Abdul-Mageed et al., 2021), (2) AraBERTv2 Antoun et al. (2020),
(3) CamelBERT-MSA and (4) CamelBERT-Mix Inoue et al. (2021). Each model is
finetuned twice for this subtask, once for predicting the corresponding
AySGNS  embedding for each input definition, and the other time for predicting
the corresponding  AyELECTRA  embedding. The final representation is computed
by taking the embedding of the CLS and passing it through a two-layer dense
network with a  AyTanh  activation function in between. The model is trained
by optimizing the Mean Squared Error (MSE) between the ground-truth
representation and the predicted one. For the learning rate scheduling policy,
we used OneCycleLR (Smith and Topin, 2017). Throughout the finetuning process,
we evaluate on the devset after every epoch, and take the checkpoint with the
highest cosine similarity score. Table 3 shows the values of the
hyperparameters used during finetuning.

To identify the optimal ensemble of our fine-tuned models, we select the model
combination that exhibited the highest performance on the devset, determined
by the cosine similarity metric, as our final solution. Tables 5 and 6 in the
Appendix shows the performance of all model combinations on the devset. The
final representation of each ensemble is taken by averaging the predicted
embedding of each for a given input definition.

Hyperparameter | Value
---|---
Batch Size | 100
lr | 1.0e-4
Learning Rate Sched. | OneCycleLR
$pct$ | 0.2
$f_{initial}$ | 25
$f_{final}$ | 100
Weight Decay | 1.0e-4
Epochs | 20
Optimizer | AdamW
Table 3: Hyperparameters Used

###  3.2 Subtask 2: English Definitions to Arabic Embeddings

Subtask 2 differs from subtask 1 by utilizing an English definition as input
instead of Arabic, with the objective of generating the embedding
representation of the Arabic word as output. Several approaches were explored
in pursuit of optimizing the system for superior output embedding quality.

#### Cross-Lingual Alignment

This method involves a two-step learning process. First, we leverage the
English Language Dictionary to learn to generate the English embeddings from
their corresponding English definition. Then the second stage utilize the
English Arabic Mapped Dictionary to learn an alignment function between both
language representations. Figure 1 shows an illustration of this model. The
motivation behind this is that the English pretrained models often yield
superior representations compared to their Arabic counterparts due to their
training on larger corpora. Here, we used RoBERTa (Liu et al., 2019) to obtain
English embeddings, following the same procedure as in subtask 1, and then
utilizing an autoencoder model to transform these embeddings into their Arabic
representations. Both the encoder and the decoder of the Autoencoder consist
of two linear layers with ReLU in between. The input and output dim is 256 and
the hidden dim is 32. However, the efficacy of converting an English
representation into an Arabic one is contingent upon the quantity of aligned
data points available in the provided resources.

#### Translate-Test

Our solution for subtask 2 that yielded the best results was inspired from
Artetxe et al. (2023). In their work, they show that machine translating a
non-English test sets into English and then running inference on a monolingual
English model can exhibit superior performance compared to using a
multilingual model, such as XLM-R Conneau et al. (2020), on the original data
zeroshot. Similarly, we use the Arabic translation of the English definitions
as input to our finetuned Arabic models. This approach enables the reuse of
models and solutions that were initially developed for subtask 1 .

##  4 Results

Subtask | Embedding | MSE | Cosine | Rank | P@1 | P@10
---|---|---|---|---|---|---
Subtask 1 | Electra | 0.152 / 0.161 | 0.645 / 0.637 | 0.242 / 0.214 | 0.031 / 0.034 | 0.099 / 0.114
SGNS | 0.030 / 0.035 | 0.605 / 0.552 | 0.254 / 0.281 | 0.445 / 0.414 | 0.597 / 0.540
Subtask 2 | Electra | 0.170 / 0.180 | 0.659 / 0.624 | 0.127 / 0.204 | 0.185 / 0.120 | 0.407 / 0.355
SGNS | 0.053 / 0.048 | 0.400 / 0.387 | 0.320 / 0.372 | 0.312 / 0.316 | 0.375 / 0.389
Table 4: Results on TestSet / DevSet for Both Subtasks. MSE is Mean-Squared-
Error. P is Precision.

Table 4 displays the results obtained on the test set across all metrics
reported in the shared task. Interestingly, the best ensemble on both subtasks
was done by taking the average of the CamelBERT-MSA and MARBERTv2 output
embeddings.

###  4.1 Subtask 1

Table 5 shows the results on the devset that we can use for further analysis.
It clearly illustrates that ensembles, regardless of the combination, enhance
the scores in comparison to using individual models. Furthermore, it is
evident that results involving CamelBERT-Mix tend to be less favorable than
those involving CamelBERT-MSA. This observation aligns with the dataset’s
nature, which predominantly features MSA definitions, thus minimizing
dialectal content.

Through examining the scores of ensembles and systems incorporating MARBERTv2
compared to those that do not, we can conclude that MARBERTv2 stands out as
the most effective model to employ or include in an ensemble among all the
tested Arabic pretrained transformers.

###  4.2 Subtask 2

The findings from Subtask 1 are applicable to Subtask 2, and this consistency
can be attributed to the reuse of models initially developed in Subtask 1 for
Subtask 2.

##  5 Discussion

#### Exploring Cross-Lingual Alignment Further

In the pursuit of optimizing our approach for the Arabic Reverse Dictionary
shared task, we implemented a cross-lingual alignment method, as detailed in
section 3.2. This method allowed us to bridge the gap between English and
Arabic definitions, by leveraging the aligned dictionary provided as part of
the shared task. Further exploration and refinement could yield promising
results in that direction.

#### Augmenting Training Data Through Self-Synthesis

In another set of experiments, we explored a very different approach that
requires further investigation in future work. The idea is to finetune of an
encoder-decoder model, such as AraT5 Nagoudi et al. (2022), jointly on two
interconnected tasks. The first task involves predicting the word embeddings
from the encoder side, while the second task entails predicting the
corresponding definition on the decoder side based on an the input word. This
approach presents an intriguing opportunity to generate diverse definition-
embedding pairs using a single model, which could subsequently be harnessed
for more robust finetuning. This self-synthesis approach could potentially
lead to better system performance by expanding the training set.

##  6 Conclusion

In this paper, we present our winning solution to the Arabic Reverse
Dictionary shared task. The objective is to derive an Arabic word
representation based on a provided definition, which can be in either Arabic
or English.

Our approach simply leverages several language models pretrained on Arabic
datasets. Through finetuning and ensembling the trained models, our method is
capable of capturing the underlying semantics of the input definitions as well
as correcting small errors done by single models.

For the first subtask, we achieve the best results by fine-tuning four Arabic
pretrained language models twice, one for predicting the  AyElectra  embedding
and once for the  AySGNS  one. This involves minimizing the discrepancy
between the predicted embedding and the model’s final representation using an
MSE loss function.

In the second subtask, our most effective solution is to repurpose the models
initially developed for the first subtask by translating the English test set
definitions into Arabic.

## Limitations

One notable limitation is related to the second subtask, where our approach
involves translating English definitions to Arabic. The results of this paper
used the existing Arabic translations that comes English test set. Therefore,
we have not investigated the quality of machine translation models, which can
significantly influence the system’s effectiveness, as inaccuracies or nuances
lost in translation may affect results. Moreover, the generalization of our
models to broader or different distributions may be constrained, as they are
optimized on specific datasets. To achieve wider applicability, we might
necessitate further finetuning on more diverse data sources. Furthermore, our
choice of evaluation metrics can influence the perceived performance of the
system, and different metrics may reveal varying aspects of its utility in
practical applications. It is essential to consider these limitations when
assessing the robustness and adaptability of our approach.

## References

  * Abdul-Mageed et al. (2021) Muhammad Abdul-Mageed, AbdelRahim Elmadany, and El Moatez Billah Nagoudi. 2021.  ARBERT & MARBERT: Deep bidirectional transformers for Arabic.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 7088–7105, Online. Association for Computational Linguistics.
  * Antoun et al. (2020) Wissam Antoun, Fady Baly, and Hazem Hajj. 2020.  Arabert: Transformer-based model for arabic language understanding.  In _LREC 2020 Workshop Language Resources and Evaluation Conference 11–16 May 2020_ , page 9.
  * Artetxe et al. (2023) Mikel Artetxe, Vedanuj Goswami, Shruti Bhosale, Angela Fan, and Luke Zettlemoyer. 2023.  Revisiting machine translation for cross-lingual classification.  _ArXiv_ , abs/2305.14240.
  * Bilac et al. (2004) Slaven Bilac, Wataru Watanabe, Taiichi Hashimoto, Takenobu Tokunaga, and Hozumi Tanaka. 2004.  Dictionary search based on the target word descrip tion. _In Proceedings of NLP_.
  * Brown and McNeill (1966) Roger Brown and David McNeill. 1966.  The “tip of the tongue” phenomenon.  _Journal of Verbal Learning and Verbal Behavior_ , 5(4):325–337.
  * Chen et al. (2018) Muhao Chen, Yingtao Tian, Haochen Chen, Kai-Wei Chang, Steven Skiena, and Carlo Zaniolo. 2018.  Learning to represent bilingual dictionaries.  _CoRR_ , abs/1808.03726.
  * Clark et al. (2020) Kevin Clark, Minh-Thang Luong, Quoc V. Le, and Christopher D. Manning. 2020.  ELECTRA: pre-training text encoders as discriminators rather than generators.  _CoRR_ , abs/2003.10555.
  * Conneau et al. (2020) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2020.  Unsupervised cross-lingual representation learning at scale.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 8440–8451, Online. Association for Computational Linguistics.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Hill et al. (2015) Felix Hill, Kyunghyun Cho, Anna Korhonen, and Yoshua Bengio. 2015.  Learning to understand phrases by embedding the dictionary.  _CoRR_ , abs/1504.00548.
  * Inoue et al. (2021) Go Inoue, Bashar Alhafni, Nurpeiis Baimukan, Houda Bouamor, and Nizar Habash. 2021\.  The interplay of variant, size, and task type in Arabic pre-trained language models.  In _Proceedings of the Sixth Arabic Natural Language Processing Workshop_ , Kyiv, Ukraine (Online). Association for Computational Linguistics.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized BERT pretraining approach.  _CoRR_ , abs/1907.11692.
  * Mikolov et al. (2013) Tomás Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean. 2013.  Efficient estimation of word representations in vector space.  In _1st International Conference on Learning Representations, ICLR 2013, Scottsdale, Arizona, USA, May 2-4, 2013, Workshop Track Proceedings_.
  * Nagoudi et al. (2022) El Moatez Billah Nagoudi, AbdelRahim Elmadany, and Muhammad Abdul-Mageed. 2022.  AraT5: Text-to-text transformers for Arabic language generation.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 628–647, Dublin, Ireland. Association for Computational Linguistics.
  * Shaw et al. (2013) Ryan Shaw, Anindya Datta, Debra VanderMeer, and Kaushik Dutta. 2013.  Building a scalable database-driven reverse dictionary.  _IEEE Transactions on Knowledge and Data Engineering_ , 25(3):528–540.
  * Smith and Topin (2017) Leslie N. Smith and Nicholay Topin. 2017.  Super-convergence: Very fast training of residual networks using large learning rates.  _CoRR_ , abs/1708.07120.
  * Yan et al. (2020) Hang Yan, Xiaonan Li, and Xipeng Qiu. 2020.  BERT for monolingual and cross-lingual reverse dictionary.  _CoRR_ , abs/2009.14790.
  * Zhang et al. (2019) Lei Zhang, Fanchao Qi, Zhiyuan Liu, Yasheng Wang, Qun Liu, and Maosong Sun. 2019\.  Multi-channel reverse dictionary model.  _CoRR_ , abs/1912.08441.

##  Appendix A Cross-Lingual Alignment Model

Figure 1: Method Explored for Subtask 2 (Section 3.2)

##  Appendix B Results on Devset

Please refer to next page.

Models | Electra | SGNS
---|---|---
MSE | Cosine | Rank | MSE | Cosine | Rank
arabert | 0.1695 | 0.6124 | 0.2491 | 0.0368 | 0.4853 | 0.3321
camelbert-mix | 0.1689 | 0.6134 | 0.2667 | 0.0389 | 0.4911 | 0.3221
camelbert-msa | 0.1681 | 0.6166 | 0.2421 | 0.0379 | 0.4947 | 0.3213
marbert | 0.1661 | 0.6265 | 0.2108 | 0.0370 | 0.5485 | 0.2782
arabert,camelbert-mix | 0.1656 | 0.6228 | 0.2525 | 0.0360 | 0.5045 | 0.3237
arabert,camelbert-msa | 0.1650 | 0.6247 | 0.2400 | 0.0358 | 0.5052 | 0.3235
arabert,marbert | 0.1618 | 0.6355 | 0.2175 | 0.0337 | 0.5511 | 0.2836
camelbert-mix,camelbert-msa | 0.1653 | 0.6239 | 0.2496 | 0.0370 | 0.5036 | 0.3208
camelbert-mix,marbert | 0.1622 | 0.6341 | 0.2267 | 0.0348 | 0.5502 | 0.2817
camelbert-msa,marbert | 0.1614 | 0.6365 | 0.2144 | 0.0345 | 0.5519 | 0.2812
arabert,camelbert-mix,camelbert-msa | 0.1642 | 0.6272 | 0.2455 | 0.0357 | 0.5095 | 0.3221
arabert,camelbert-mix,marbert | 0.1616 | 0.6356 | 0.2286 | 0.0339 | 0.5466 | 0.2862
arabert,camelbert-msa,marbert | 0.1610 | 0.6371 | 0.2204 | 0.0338 | 0.5472 | 0.2860
camelbert-mix,camelbert-msa,marbert | 0.1614 | 0.6361 | 0.2268 | 0.0346 | 0.5452 | 0.2849
arabert,camelbert-mix,camelbert-msa,marbert | 0.1613 | 0.6363 | 0.2287 | 0.0341 | 0.5421 | 0.2895
Table 5: Performance Analysis on the Devset of Subtask-1 Using Various Model Ensembles. Models | Electra | SGNS
---|---|---
MSE | Cosine | Rank | MSE | Cosine | Rank
arabert | 0.1879 | 0.6014 | 0.2369 | 0.0491 | 0.3500 | 0.3925
camelbert-mix | 0.1894 | 0.5974 | 0.2482 | 0.0520 | 0.3504 | 0.3956
camelbert-msa | 0.1860 | 0.6066 | 0.2167 | 0.0498 | 0.3580 | 0.3845
marbert | 0.1858 | 0.6108 | 0.2115 | 0.0530 | 0.3818 | 0.3739
arabert,camelbert-mix | 0.1848 | 0.6104 | 0.2354 | 0.0486 | 0.3619 | 0.3920
arabert,camelbert-msa | 0.1829 | 0.6149 | 0.2218 | 0.0479 | 0.3640 | 0.3855
arabert,marbert | 0.1806 | 0.6226 | 0.2106 | 0.0478 | 0.3867 | 0.3752
camelbert-mix,camelbert-msa | 0.1842 | 0.6117 | 0.2245 | 0.0494 | 0.3619 | 0.3894
camelbert-mix,marbert | 0.1821 | 0.6186 | 0.2175 | 0.0493 | 0.3838 | 0.3776
camelbert-msa,marbert | 0.1800 | 0.6238 | 0.2038 | 0.0484 | 0.3874 | 0.3715
arabert,camelbert-mix,camelbert-msa | 0.1827 | 0.6160 | 0.2248 | 0.0481 | 0.3662 | 0.3890
arabert,camelbert-mix,marbert | 0.1808 | 0.6222 | 0.2162 | 0.0476 | 0.3841 | 0.3778
arabert,camelbert-msa,marbert | 0.1794 | 0.6255 | 0.2075 | 0.0472 | 0.3860 | 0.3736
camelbert-mix,camelbert-msa,marbert | 0.1805 | 0.6228 | 0.2135 | 0.0482 | 0.3836 | 0.3761
arabert,camelbert-mix,camelbert-msa,marbert | 0.1800 | 0.6240 | 0.2125 | 0.0474 | 0.3830 | 0.3782
Table 6: Performance Analysis on the Devset of Subtask-2 Using Various Model
Ensembles.


# Fitting ellipses to noisy measurements

Sebastian Dingler1 (orcid: 0000-0002-0162-8428) 1This paper is an independent
contribution without affiliation. To mail me remove foobar
<EMAIL_ADDRESS>
###### Abstract

This work deals with fitting of ellipses to noisy measurements. The literature
knows many different approaches for this. The main representatives are
presented and discussed in this paper. Furthermore, the case is considered
when outliers are present in the measurement data. Robust methods which are
less sensitive to outliers are suitable for this case. All discussed methods
are compared by a simulation. The code for the simulation is available for
free use on github.com/sebdi/ellipse-fitting.

##  I Introduction

Approximating given data using a model is a typical task in the context of
pattern recognition, tracking and image processing. The goal is to reduce the
complexity of the data and to provide higher application layers with a simpler
representation of the data. Occasionally it is essential to estimate a simple
model with measurements. As an example the tracking of objects can be named.
In case of high measurement noise, few measurements or objects with a strong
motion it is difficult to estimate the shape of an object. In this case it is
necessary to use a simple shape with few degrees of freedom to describe the
object. For this purpose the ellipse is one of the most used models. In image
processing ellipses often appear because a circle becomes an ellipse by
perspective projection [1, 2, 3]. The same is true for spherical objects [4].
Moreover, because of the non-linearity of the ellipse, it is a suitable model
problem to test and compare different regression algorithms. Due to the non-
linearity, several methods exist to fit ellipses to measurements. Each method
uses a different assumption or simplification. In this paper, these methods
will be presented.

To this end, Section II lays the necessary groundwork. Sections III through VI
then present the various methods. A simulation in Section VII and a summary
with outlook in Section VIII complete this work.

##  II Mathematical description of ellipses

An ellipse $C$ can be described for the coordinate axes $x$ and $y$ as follows

$C(x,y)={\mathbf{a}}x^{2}+2{\mathbf{b}}xy+{\mathbf{c}}y^{2}+2{\mathbf{d}}x+2{\mathbf{e}}y+{\mathbf{f}}=0\enspace.$
(1)

Here ${\mathbf{a}}$, ${\mathbf{b}}$ and ${\mathbf{c}}$ cannot be zero at the
same time. Another constraint is that
${\mathbf{b}}^{2}-{\mathbf{a}}{\mathbf{c}}<0$ must hold. Given noisy
measurements for this ellipse, the goal is to determine the parameters
${\mathbf{a}}$, ${\mathbf{b}}$, ${\mathbf{c}}$, ${\mathbf{d}}$, ${\mathbf{e}}$
and ${\mathbf{f}}$ in such a way that they best describe the data in terms of
a quality criterion. Because the trivial solution
${\mathbf{a}}={\mathbf{b}}={\mathbf{c}}={\mathbf{d}}={\mathbf{e}}={\mathbf{f}}=0$
must be prevented, the literature often uses three normalizations. These are:

  * •

Normalization with ${\mathbf{a}}+{\mathbf{c}}=1$

  * •

Normalization with
${\mathbf{a}}^{2}+{\mathbf{b}}^{2}+{\mathbf{c}}^{2}+{\mathbf{d}}^{2}+{\mathbf{e}}^{2}+{\mathbf{f}}^{2}=1$

  * •

Normalization with ${\mathbf{f}}=1$.

Another useful form of the elliptic equation is the general parametric form

$\displaystyle x(t)$ $\displaystyle=x_{c}+a\ \text{cos}(t)\
\text{cos}(\alpha)-b\ \text{sin}(t)\ \text{sin}(\alpha)\enspace,$ (2)
$\displaystyle y(t)$ $\displaystyle=y_{c}+a\ \text{cos}(t)\
\text{sin}(\alpha)+b\ \text{sin}(t)\ \text{cos}(\alpha)\enspace,$ (3)

where the parameter $t\in[0,2\pi]$ can be varied to get different values of
the ellipse. The other parameters are the center of the ellipse with $x_{c}$
and $y_{c}$, the semi-major axis $a$, the semi-minor axis $b$ and $\alpha$ the
angle between $x$-axis and the semi-major axis $a$.

Furthermore, an ellipse can also be represented without rotation and a center
at the origin of the coordinate system. This coordinate system is denoted with
$XY$ to distinguish it from the general coordinate system $xy$. The
transformation between the coordinate systems $xy$ and $XY$ is calculated with
the rotation matrix

$\displaystyle{\mathbf{R}}=\begin{bmatrix}\text{cos}\alpha&\text{sin}\alpha\\\\[3.00003pt]
-\text{sin}\alpha&\text{cos}\alpha\end{bmatrix}\enspace.$ (4)

The transformation for a point ${\mathbf{x}}=[x,y]^{T}$ in the $xy$ coordinate
system is then done with

$\displaystyle{\mathbf{X}}={\mathbf{R}}({\mathbf{x}}-{\mathbf{x}}_{c})\enspace,$
(5)

with ${\mathbf{x}}_{c}=[x_{c},y_{c}]^{T}$ and ${\mathbf{X}}=[X,Y]^{T}$. The
reverse transformation with

$\displaystyle{\mathbf{x}}={\mathbf{R}}^{-1}{\mathbf{X}}+{\mathbf{x}}_{c}\enspace.$
(6)

##  III Least-squares fitting based on algebraic distance

In this section, the least-squares method will be used to show how the
parameters of the ellipse can be determined based on noisy measurements.
Applying the normalization ${\mathbf{a}}+{\mathbf{c}}=1$ to the elliptic
equation, this eliminates ${\mathbf{c}}$ and leads to

$\displaystyle
C(x,y)={\mathbf{a}}(x^{2}-y^{2})+2{\mathbf{b}}xy+2{\mathbf{d}}x+2{\mathbf{e}}y+y^{2}+{\mathbf{f}}\enspace.$
(7)

This can be converted into vector notation with
${\mathbf{\phi}}=[x^{2}-y^{2},2xy,2x,2y,1]$,
${\mathbf{p}}=[{\mathbf{a}},{\mathbf{b}},{\mathbf{d}},{\mathbf{e}},{\mathbf{f}}]^{T}$
and $z=-y^{2}$, to

$\displaystyle C(x,y)={\mathbf{\phi}}{\mathbf{p}}-z.$ (8)

If $n$ measurements
$\left\\{(x_{1},y_{1}),(x_{2},y_{2}),....,(x_{n},y_{n})\right\\}$ are given,
this leads to $n$ equations

$\displaystyle C_{1}(x_{1},y_{1})$ $\displaystyle=$
$\displaystyle{\mathbf{\phi}}_{1}{\mathbf{p}}-z_{1}$ (9) $\displaystyle
C_{2}(x_{2},y_{2})$ $\displaystyle=$
$\displaystyle{\mathbf{\phi}}_{2}{\mathbf{p}}-z_{2}$ (10)
$\displaystyle\vdots$ $\displaystyle=$ $\displaystyle\vdots$ (11)
$\displaystyle C_{n}(x_{n},y_{n})$ $\displaystyle=$
$\displaystyle{\mathbf{\phi}}_{n}{\mathbf{p}}-z_{n}\enspace.$ (12)

These can be expressed in vector notation. For this purpose the following is
defined

$\displaystyle{\mathbf{C}}$ $\displaystyle=$
$\displaystyle[C_{1}(x_{1},y_{1}),C_{2}(x_{2},y_{2}),..,C_{n}(x,y)]^{T}\enspace,$
(13) $\displaystyle{\mathbf{\Phi}}$ $\displaystyle=$
$\displaystyle[{\mathbf{\phi}}_{1}^{T},{\mathbf{\phi}}_{2}^{T},..,{\mathbf{\phi}}_{n}^{T}]^{T}\enspace,$
(14) $\displaystyle{\mathbf{y}}$ $\displaystyle=$
$\displaystyle[z_{1},z_{2},...,z_{n}]^{T}\enspace.$ (15)

This yields the equation

${\mathbf{C}}={\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}}\enspace.$ (16)

The goal is now to find a ${\mathbf{p}}$ for noisy measurements that optimally
represents these measurements with respect to the least-squares criterion.
Since measurements are always noisy, there will usually be no ${\mathbf{p}}$
that is zero for all measurements ${\mathbf{C}}$, hence an error
$\tilde{{\mathbf{e}}}$ with
$\tilde{{\mathbf{e}}}=[\tilde{e}_{1},\tilde{e}_{2},...,\tilde{e}_{n}]^{T}$
will remain. Therefore, the real goal is to minimize the error
$\tilde{{\mathbf{e}}}$ by choosing ${\mathbf{p}}$ appropriately. For this
purpose, linear regression finds a vector ${\mathbf{p}}$ by minimizing the
quality function

$\displaystyle{\mathbf{G}}({\mathbf{p}})$ $\displaystyle=$
$\displaystyle{\mathbf{C}}^{T}{\mathbf{C}}$ (17) $\displaystyle=$
$\displaystyle\tilde{{\mathbf{e}}}^{T}\tilde{{\mathbf{e}}}$ (18)
$\displaystyle=$
$\displaystyle({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})^{T}({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})\enspace.$
(19)

This is achieved by setting the derivative of ${\mathbf{G}}$ equal to zero,
resulting in

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial{\mathbf{p}}}{\mathbf{G}}({\mathbf{p}})$
(20) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial{\mathbf{p}}}({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})^{T}({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})$
(21) $\displaystyle 0$ $\displaystyle=$ $\displaystyle
2{\mathbf{\Phi}}^{T}({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})$ (22)
$\displaystyle{\mathbf{p}}$ $\displaystyle=$
$\displaystyle({\mathbf{\Phi}}^{T}{\mathbf{\Phi}})^{-1}{\mathbf{\Phi}}^{T}{\mathbf{y}}\enspace.$
(23)

Thus ${\mathbf{p}}$ is the solution with the smallest squared error with
respect to the algebraic distance since $C(x,y)$ is minimized.

##  IV Least-squares fitting based on orthogonal distance

In Section III the algebraic distance $C(x,y)$ was minimized. This resulted in
a very lean and easy to compute solution for ${\mathbf{p}}$. However, the
least-squares solution based on the algebraic distance is not entirely
correct. The error made in this method is called the high curvature bias.
Explained in simple terms, the bias arises from the given shape of the
ellipse. Thus, a measurement along the semi-major axis $a$ has a large value
for $C(x,y)$ while a measurement along the semi-minor axis $b$ has a smaller
value. However, both measurements can be equally close to the ellipse although
their distance $C(x,y)$ is different. To avoid this bias, the orthogonal
distance between a measurement $(x_{i},y_{i})$ and the ellipse must be used
(cf. Figure 1). The orthogonal distance $d$ is the shortest connection between
a point $(x_{i},y_{i})$ and the ellipse. The parameter vector ${\mathbf{p}}$
can be obtained by minimizing the sum of the squared distances

${\mathbf{G}}({\mathbf{p}})=\sum^{n}_{i=1}d^{2}_{i}\enspace.$ (24)

For the 2-dimensional case, the distance $d$ is a 4th degree polynomial and
can be solved analytically [5, 6]. However, according to [3] the analytic
solution is unstable for measurements $|x_{i}|\approx 0$ or $|y_{i}|\approx
0$, so in the following section the solution of [3] which uses the Gauss-
Newton algorithm will be presented. Once the orthogonal distance $d_{i}$ for
each point $i$ is found, the quality function ${\mathbf{G}}$ must be minimized
to find the parameter vector ${\mathbf{p}}$ with the smallest error. For this
purpose the iterative algorithms like: Gauss-Newton, Steepest Gradient
Descent, Levenberg-Marquardt or the Simplex algorithm are appropriate. In
Section IV-B, the Gauss-Newton approach from [3] will be used.

$(x_{i},y_{i})$$(x_{c},y_{c})$$d_{i}$$(x_{\bot},y_{\bot})$(0,1.22)(0.8,0.8)-1175
$x$$y$ Figure 1: Orthogonal distance $d_{i}$ between a measurement
$(x_{i},y_{i})$ and the orthogonal point $(x_{\bot},y_{\bot})$ on the ellipse.

###  IV-A Determination of the orthogonal point on the ellipse

To get the orthogonal distance $d$ the corresponding orthogonal point
$(x_{\bot},y_{\bot})$ on the ellipse is calculated. The calculation is
simplified if this is done in the $XY$ coordinate system, i.e. in the
coordinate system in which the ellipse has no rotation and is symmetrical to
each coordinate axis. The elliptic equation for this case is

$\displaystyle\frac{X^{2}}{a^{2}}+\frac{Y^{2}}{b^{2}}=1$ (25)

or transformed

$\displaystyle f_{1}(X,Y)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(a^{2}Y^{2}+b^{2}X^{2}-a^{2}b^{2})=0\enspace.$ (26)

For the condition that the tangent of the orthogonal point
$(X_{\bot},Y_{\bot})$ is orthogonal to a measurement $(X_{i},Y_{i})$

$\displaystyle\frac{dy}{dx}\cdot\frac{Y_{i}-Y_{\bot}}{X_{i}-X_{\bot}}=\frac{-b^{2}X_{\bot}}{a^{2}Y_{\bot}}\cdot\frac{Y_{i}-Y_{\bot}}{X_{i}-X_{\bot}}=-1$
(27)

the second equation

$\displaystyle f_{2}(X,Y)$ $\displaystyle=$ $\displaystyle
b^{2}X(Y_{i}-Y)-a^{2}Y(X_{i}-X)=0$ (28)

can be derived. These two equations

$\displaystyle{\mathbf{f}}(X,Y)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}f_{1}(X,Y)\\\\[3.00003pt]
f_{2}(X,Y)\end{bmatrix}$ (29)

have to be solved in vector notation. Although this is analytically possible,
it exhibits numerical instability as already described. Therefore, in the
following the Gauss-Newton method shall be applied to the system of equations.
For this the Jacobian matrix

$\displaystyle{\mathbf{Q}}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\frac{\partial f_{1}}{\partial x}&\frac{\partial
f_{1}}{\partial y}\\\\[3.00003pt] \frac{\partial f_{2}}{\partial
x}&\frac{\partial f_{2}}{\partial y}\end{bmatrix}$ (30) $\displaystyle=$
$\displaystyle\begin{bmatrix}b^{2}x&a^{2}y\\\\[3.00003pt]
(a^{2}-b^{2})y+b^{2}y_{i}&(a^{2}-b^{2})x-a^{2}x_{i}\end{bmatrix}$ (31)

is calculated. The initial value ${\mathbf{X}}_{0}$ can be calculated by
averaging the points ${\mathbf{X}}_{k1}$ and ${\mathbf{X}}_{k2}$.

$\displaystyle{\mathbf{X}}_{0}=\frac{1}{2}({\mathbf{X}}_{k1}+{\mathbf{X}}_{k2})$
(32)

with

$\displaystyle{\mathbf{X}}_{k1}=\frac{\begin{bmatrix}X_{i}\\\\[3.00003pt]
Y_{i}\end{bmatrix}ab}{\sqrt{b^{2}X_{i}^{2}+a^{2}Y^{2}_{i}}}$ (33)

and

$\displaystyle{\mathbf{X}}_{k2}=\begin{cases}\begin{bmatrix}X_{i}\\\\[3.00003pt]
\text{sign}(Y_{i})\frac{b}{a}\sqrt{a^{2}-X_{i}^{2}}\end{bmatrix}&\text{if}\
|X_{i}|<a,\\\ \begin{bmatrix}\text{sign}(X_{i})a\\\\[3.00003pt]
0\end{bmatrix}&\text{if}\ |X_{i}|\geq a\enspace.\end{cases}$ (34)

The orthogonal point ${\mathbf{X}}^{\prime}=[X_{\bot},Y_{\bot}]^{T}$ is
calculated iterating

$\displaystyle{\mathbf{X}}^{\prime}_{k+1}={\mathbf{X}}^{\prime}_{k}-{\mathbf{Q}}^{-1}_{k}{\mathbf{f}}(X_{k})\enspace.$
(35)

Finally, by transforming back into the $xy$ coordinate system, the sought
vector of the orthogonal distance can be determined

$\displaystyle{\mathbf{d}}_{i}={\mathbf{x}}_{i}-{\mathbf{x}}^{\prime}_{i}\enspace.$
(36)

###  IV-B Fitting an ellipse with orthogonal distances

With the orthogonal point ${\mathbf{X}}^{\prime}$ from section IV-A it is now
possible to perform the fitting of the ellipse also with the Gauss-Newton
algorithm. For this purpose the vector
${\mathbf{q}}=[x_{c},y_{c},a,b,\alpha]^{T}$ is defined which contains the
searched parameters. For each measurement ${\mathbf{x}}_{i}$ a Jacobi matrix

$\displaystyle{\mathbf{J}}_{{\mathbf{x}}={\mathbf{x}}^{\prime}_{i}}=({\mathbf{R}}^{-1}{\mathbf{Q}}^{-1}{\mathbf{B}})$
(37)

gets computed with

$\displaystyle{\mathbf{B}}$ $\displaystyle=$ $\displaystyle\begin{bmatrix}\
{\mathbf{B}}_{1}\ {\mathbf{B}}_{2}\ {\mathbf{B}}_{3}\
{\mathbf{B}}_{4}\end{bmatrix}\enspace,$ (38) $\displaystyle{\mathbf{B}}_{1}$
$\displaystyle=$ $\displaystyle\begin{bmatrix}b^{2}xC-a^{2}yS\\\\[3.00003pt]
b^{2}(yi-y)C+a^{2}(xi-x)S\end{bmatrix}\enspace,$ (39)
$\displaystyle{\mathbf{B}}_{2}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}b^{2}xS+a^{2}yC\\\\[3.00003pt]
b^{2}(yi-y)S-a^{2}(xi-x)C\end{bmatrix}\enspace,$ (40)
$\displaystyle{\mathbf{B}}_{3}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}a(b^{2}-y^{2})\\\\[3.00003pt]
2ay(xi-x)\end{bmatrix}\enspace,$ (41) $\displaystyle{\mathbf{B}}_{4}$
$\displaystyle=$ $\displaystyle\begin{bmatrix}b(a^{2}-x^{2})\\\\[3.00003pt]
-2bx(yi-y)\end{bmatrix}\enspace.$ (42)

The sought parameter vector ${\mathbf{q}}$ is then iteratively determined with
a step size $\lambda$ using

$\displaystyle{\mathbf{q}}_{k+1}={\mathbf{q}}_{k}+\lambda\Delta{\mathbf{q}}$
(44)

and

$\displaystyle\begin{bmatrix}{\mathbf{J}}_{{\mathbf{x}}={\mathbf{x}}^{\prime}_{1}}\\\\[3.00003pt]
{\mathbf{J}}_{{\mathbf{x}}={\mathbf{x}}^{\prime}_{2}}\\\\[3.00003pt]
\vdots\\\\[3.00003pt]
{\mathbf{J}}_{{\mathbf{x}}={\mathbf{x}}^{\prime}_{n}}\end{bmatrix}\Delta{\mathbf{q}}=\begin{bmatrix}{\mathbf{d}}_{1}\\\\[3.00003pt]
{\mathbf{d}}_{2}\\\\[3.00003pt] \vdots\\\\[3.00003pt]
{\mathbf{d}}_{n}\end{bmatrix}.$ (45)

##  V Gradient weighted least-squares fitting

For the fitting of ellipses the least-squares solution based on the algebraic
distance from Section III provides satisfying results. If accuracy is
required, the least-squares solution based on the orthogonal distance from
Section IV-A should be chosen. However, the high accuracy is bought with a
considerably higher computational effort. Thus, for each measurement, the
orthogonal distance must first be found before iteratively determining the
solution for ${\mathbf{q}}$. In this section, the Gradient Weighted-Least-
Squares method will be presented, which is much faster than least-squares
based on orthogonal distance and yet more accurate than the least-squares
method based on algebraic distance. The Gradient Weighted-Least-Squares uses a
normalization with the help of the gradient which leads to a so-called
statistically optimal weighted algebraic fit [5]. It is sufficient to divide
the original function with its gradient. Thus, similar to Section III, the
quadratic goodness measure

$\displaystyle{\mathbf{G}}({\mathbf{p}})$ $\displaystyle=$ $\displaystyle\sum
C^{{}^{\prime}2}_{i}(x_{i},y_{i})$ (46) $\displaystyle=$
$\displaystyle\sum\frac{C^{2}_{i}(x_{i},y_{i})}{||\nabla
C_{i}(x_{i},y_{i})||^{2}}$ (47) $\displaystyle=$
$\displaystyle({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})^{T}{\mathbf{W}}^{-1}({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})$
(48)

can be minimized. Where the normalization is done in the weighting matrix
${\mathbf{W}}$, which is

${\mathbf{W}}=\text{diag}(||\nabla C_{1}||^{2},||\nabla
C_{2}||^{2},...,||\nabla C_{n}||^{2})\enspace.$ (49)

The gradient is calculated from

$||\nabla C_{i}(x_{i},y_{i})||^{2}=\bigg{(}\frac{\partial C_{i}}{\partial
x_{i}}\bigg{)}^{2}+\bigg{(}\frac{\partial C_{i}}{\partial y_{i}}\bigg{)}^{2}$
(50)

with

$\displaystyle\frac{\partial C_{i}}{\partial x_{i}}$ $\displaystyle=$
$\displaystyle 2{\mathbf{a}}x_{i}+2{\mathbf{b}}y_{i}+2{\mathbf{d}}\enspace,$
(51) $\displaystyle\frac{\partial C_{i}}{\partial y_{i}}$ $\displaystyle=$
$\displaystyle-2{\mathbf{a}}y_{i}+2{\mathbf{b}}x_{i}+2{\mathbf{e}}+2{\mathbf{y}}\enspace.$
(52)

This is known as the weighted least-squares method which is calculated by
zeroing the first derivative to

$\displaystyle 0$
$\displaystyle=\frac{\partial}{\partial{\mathbf{p}}}{\mathbf{G}}({\mathbf{p}})$
(53) $\displaystyle 0$
$\displaystyle=\frac{\partial}{\partial{\mathbf{p}}}({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})^{T}{\mathbf{W}}^{-1}({\mathbf{X}}{\mathbf{p}}-{\mathbf{y}})$
(54) $\displaystyle 0$
$\displaystyle=2{\mathbf{\Phi}}^{T}{\mathbf{W}}^{-1}({\mathbf{\Phi}}{\mathbf{p}}-{\mathbf{y}})$
(55) $\displaystyle{\mathbf{p}}$
$\displaystyle=({\mathbf{\Phi}}^{T}{\mathbf{W}}^{-1}{\mathbf{\Phi}})^{-1}{\mathbf{\Phi}}^{T}{\mathbf{W}}^{-1}{\mathbf{y}}\enspace.$
(56)

Thus the solution for ${\mathbf{p}}$ looks at first sight like a closed
solution for ${\mathbf{p}}$. Unfortunately, this is not the case since
${\mathbf{W}}$ depends on ${\mathbf{p}}$. Therefore, an iterative procedure
must be used here, which is structured as follows:

  1. 1.

Calculate an initial value for ${\mathbf{p}}^{(k=0)}$ (for example, using the
least-squares solution based on the algebraic distance).

  2. 2.

Calculate the weighting matrix ${\mathbf{W}}$ based on the current value of
${\mathbf{p}}^{(k)}$ and all measurements.

  3. 3.

Compute a new ${\mathbf{p}}^{(k+1)}$ using the equation (56)

  4. 4.

If ${\mathbf{p}}^{(k)}$ and ${\mathbf{p}}^{(k+1)}$ are very similar, then
finish the procedure. If not, then continue with step 2.

##  VI Robust Estimation

In the previous sections, methods were presented which allow the fitting of
ellipses. In practice, however, all the methods mentioned have in common that
they react very sensitively to outliers in the measurement data. Outliers are
incorrect data points or measurements that occur due to various causes. In
image processing systems, for example, these are optical reflections that lead
to outliers in the feature space. Since these outliers occur unexpectedly in
practice, they cannot be modeled exactly or excluded in the conceptual design
of a system. Therefore, there is no general approach, but it must be
considered from application to application how the estimator can be stabilized
against outliers. A number of techniques exist in the literature for this
purpose, for example:

  * •

Clustering with the Hough transform: Here, the measurements are transformed
into Hough space and then a maximum likelihood estimator is used to determine
the maximum of the cluster. Since it is rarely applicable to up to more than 3
unknowns, it is therefore not suitable for fitting an ellipse.

  * •

Regression Diagnostic: Here, an initial estimate is made using an arbitrary
procedure. Then, individual measurements are removed based on the residuals if
they are too large with respect to a threshold.

In the following sections, the M-Estimator and the Least Median of Squares
will be considered, since these methods are suitable for fitting ellipses.

###  VI-A M-Estimator

The method of least-squares tries to minimize the squared error. If the
measured data have outliers, these are weighted very strongly by the quadratic
term. A single outlier can distort the result of an estimator so much that it
takes many correct measurements to compensate for this error. M-estimators
therefore use an error function $\rho$ that does not grow as much as the
quadratic error criterion. However, this property should be the only change to
the error function. The error function should still be symmetric, positive
definite, and have a minimum at zero. These requirements are met by many
functions and it is ultimately a design parameter that can be chosen for any
application (a review of these functions can be found in [6]). In general,
instead of the sum of squared distances $d$, a function $\rho$ is minimized as
a function of the distance $d$, formally expressed this is

$\text{min}\sum_{i}\rho(d_{i})\enspace.$ (57)

However, this quality function is not solved directly but via the so-called
Reweighted Least-Squares approach. In this approach a vector ${\mathbf{p}}$
with ${\mathbf{p}}=[p_{1},p_{2},...,p_{n}]^{T}$ is searched to solve the
problem of (57). The solution is obtained by setting the first derivative to
zero. Furthermore, an Influence Function with $\psi(x)=d\rho(x)/dx$ and an
Weight Function with $w(x)=\psi(x)/x$ are defined. The approach thus has $m$
equations of the form

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\sum^{m}_{i}\psi(d_{i})\frac{\partial d_{i}}{\partial p_{j}}$
(58) $\displaystyle=$ $\displaystyle\sum^{m}_{i}w(d_{i})d_{i}\frac{\partial
d_{i}}{\partial p_{j}}\enspace,$ (59)

that solves the following Iterated reweighted least-squares problem

$\displaystyle\text{min}\sum^{m}_{i}w(d^{k-1}_{i})d^{2}_{i}.$ (60)

This is done by recalculating the weight $w(d^{k-1}_{i})$ in each iteration
step $k$ based on the previous distance $d^{k-1}_{i}$ of step $k-1$. For many
functions for $\psi(x)$ and thus for $w(x)$ no closed-form solution exists and
therefore the use of an iterative algorithm is also necessary here (e.g.
Newton’s method). Furthermore, in most cases the weighted least-squares
approach can be used by performing a reweighting in each step.

TABLE I: The parameters of the estimated ellipses compared to the ideal ellipse for the simulation without outliers. method | semi-major axis $a$ | semi-minor axis $b$ | center $(x_{c},y_{c})$ | rotation $\alpha$
---|---|---|---|---
Ground truth ellipse | 24.00 | 12.00 | $(0.00,0.00)$ | 0.00
Least-squares fitting using algebraic distance | 21.98 | 11.80 | $(-0.91,0.02)$ | 3.14
Least-squares fitting using orthogonal distance | 23.33 | 11.94 | $(-0.53,-0.53)$ | 3.14
Gradient weighted least-squares fitting | 22.45 | 11.78 | $(-0.71,0.00)$ | 3.14

###  VI-B Least Median of Squares

Another way to handle outliers is to use the median and minimize it:

$\displaystyle\text{min}\ \text{Median}\ d^{2}_{i}\enspace.$ (61)

The median is known to be particularly robust against outliers. The problem is
that it is very difficult to formulate an analytical formula for the median,
since the values have to be ordered by size in order to use the median. Since
this is relatively time-consuming for a large data set, only a randomly
selected subset of the measurements is used for the estimation. Such a
procedure can look as follows:

  1. 1.

Choose $m$ random subsets $J_{i}=\left\\{(x_{j},y_{j})\right\\}^{p}_{j=1}$ of
measurements. Each subset $J_{i}$ should contain at least $p=5$ measurements,
since 5 measurements can be used to determine an ellipse.

  2. 2.

For each subset $J_{i}$, an ellipse $C_{i}$ is determined using any method.
The ellipse $C_{i}$ is described by a parameter vector ${\mathbf{p}}_{i}$.

  3. 3.

For each ellipse $C_{i}$, the median $M_{i}$ of the residuals $r$ to each
measurement of the entire data set of size $n$ is determined. Formally, this
is
$M_{i}=\underset{j=1,...,n}{\text{median}}(r^{2}_{i}({\mathbf{p}}_{i},{\mathbf{x}}_{j}))$.
A choice can be made whether to use the algebraic distance or the more
computationally expensive orthogonal distance.

  4. 4.

The ellipse $C_{i}$ and the parameter vector ${\mathbf{p}}_{i}$ with the
minimum $M_{i}$ are chosen as the best estimate.

This procedure is similar to the procedure used by the RANSAC algorithm. It
should be noted that the RANSAC algorithm is very popular and many
implementations are available. Therefore, for practical use, RANSAC can be
preferred to the Least Median of Squares.

gt.csv noise.csv LSA.csv LSO.csv LSWA.csv
[plotstyle=line,linecolor=black,linewidth=2pt,linestyle=solid,dotsep=0.4pt,yStart=-25,yEnd=25]
[plotstyle=dots,linecolor=black,linewidth=0.1pt,linestyle=dotted,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1]
[plotstyle=line,linecolor=blue,linewidth=1pt,linestyle=dashed,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1]
[plotstyle=line,linecolor=green,linewidth=1pt,linestyle=dashed,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1]
[plotstyle=line,linecolor=red,linewidth=1pt,linestyle=dotted,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1]
Gradient Weighted FittingLeast-Squares fitting using orthogonal distanceLeast-
Squares fitting using algebraic distanceNoisy MeasurementsGround Truth Ellipse
Figure 2: Qualitative result of the simulation. It can be clearly seen that
the least-squares method based on the orthogonal distance estimates the
ellipse best. The least-squares estimation based on the algebraic distance has
the highest error.

##  VII Simulation

This section aims to demonstrate the performance of the algorithms shown. The
simulation code is publicly available at github.com/sebdi/ellipse-fitting,
which allows any reader to reproduce the results or run simulations with
different parameters.

###  VII-A Fitting with noisy measurements without outliers

Noisy measurements were generated using equation (2) and (3). For this, the
parameter $t$ is varied in the interval $[\frac{\pi}{2},\frac{3\pi}{2}]$ and
homogeneous-isotropic noise with constant variance $\sigma^{2}=0.25$ is added
to $x(t)$ and $y(t)$. Then, an estimate of the ellipse is made using these
measurements. For this the presented methods

  * •

least-squares with algebraic distance,

  * •

least-squares with orthogonal distance,

  * •

and gradient weighted least-squares.

is used. Figure 2 shows the qualitative result of such a simulation. Table I
further shows the parameters of the ellipse and the estimates of each method.
The simulation shows that the best result is obtained with the least-squares
method based on the orthogonal distance. The estimate based on the algebraic
distance has the largest error, followed by the Gradient Weighted Least-
Squares method which is a good compromise between the other two methods.

###  VII-B Fitting with noisy measurements with outlier

TABLE II: The parameters of the estimated ellipses compared to the ideal ellipse for the simulation with outliers. method | semi-major axis $a$ | semi-minor axis $b$ | center $(x_{c},y_{c})$ | rotation $\alpha$
---|---|---|---|---
Ground truth ellipse | 24.00 | 12.00 | (0.00,0.00) | 0.00
Least-squares fitting using algebraic distance | 8.44 | 7.40 | (-5.28,-0.00) | 0.00
Least-squares fitting using orthogonal distance | 28.67 | 12.96 | (4.65,4,65) | 3.16
Gradient weighted least-squares fitting | 6.94 | 6.63 | (-5.39,0.24) | 0.00
M-Estimator (Cauchy) | 21.56 | 11.76 | (-1.12,0.01) | 0.00
Least median of squares | 22.56 | 11.71 | (-0.625,0.13) | 3.13

As mentioned earlier, measured data often have outliers that are not accounted
for in the standard procedures from Section III to V. In this simulation, some
outliers were added to the measurements to highlight the problem. In addition
to the previous methods, the methods

  * •

M-Estimator

  * •

Least Median of Squares

were simulated, that take outliers into account. The qualitative result is
shown in Fig. 3 to be examined. The estimated parameters of the ellipses can
be found in Tab. II. The simulation shows that the methods that take outliers
into account are much more robust compared to the standard methods.

out_gt.csv out_noise.csv out_LSA.csv out_LSO.csv out_GWLS.csv out_MEst.csv
out_LSoS.csv
[plotstyle=line,linecolor=black,linewidth=2pt,linestyle=solid,dotsep=0.4pt,yStart=-25,yEnd=25][plotstyle=dots,linecolor=black,linewidth=0.1pt,linestyle=dotted,dotsep=0.4pt,yStart=-15,yEnd=25,opacity=1][plotstyle=line,linecolor=red,linewidth=1pt,linestyle=dotted,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1][plotstyle=line,linecolor=red,linewidth=1pt,linestyle=dashed,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1][plotstyle=line,linecolor=blue,linewidth=1pt,linestyle=dashed,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1][plotstyle=line,linecolor=blue,linewidth=1pt,linestyle=dotted,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1][plotstyle=line,linecolor=green,linewidth=1pt,linestyle=dotted,dotsep=0.4pt,yStart=-25,yEnd=25,opacity=1]
Least Median of SquaresM-Estimator (Cauchy)Gradient Weighted Least Squares
FittingLeast-Squares Fitting using orthogonal distanceLeast-Squares Fitting
using algebraic distanceNoisy Measurements with OutlierGround Truth Ellipse
Figure 3: Qualitative result of the simulation with measurements that have
outliers. It can be clearly seen that methods like M-Estimator and Least
Median of Squares achieve a better result because they take into account
outliers in the measurement data.

##  VIII CONCLUSION

In this work, methods were discussed with which ellipses can be approximated
on the basis of noisy measurements. For this purpose, a very simple method was
presented with the least-square solution based on the algebraic distance.
However, this method is inaccurate because it has the so-called High Curvature
Bias. This problem is addressed by the orthogonal distance, but its solution
is more computationally expensive. Furthermore, the Gradient Weighted-Least-
Squares approach was presented as a method that, although it uses the
algebraic distance, achieves good accuracy. However, the Gradient Weighted-
Least-Squares solution is much faster to compute than the least-squares method
based on the orthogonal distance.

Furthermore, robust methods were presented that are less sensitive to outliers
in the measurements. For this purpose, the M-estimators use a different error
function that weights outliers less than the squared error measure. The other
method is to use the median where large outliers are ignored. Last but not
least, the theoretical explanations were demonstrated with the help of a
simulation.

All methods shown have in common that they consider the problem of noisy
measurements as a deterministic problem. However, there are also stochastic
methods such as [7, 8] that solve the problem using a Bayesian approach. In
practice, therefore, it must be weighed whether to use a simple and fast
method such as the least-squares approach based on algebraic distance or a
more computationally expensive approach. For this purpose, the simulation
provided on github.com can be used to make initial attempts to determine which
method is more suitable for the application problem.

Also, the methods shown do not have a recursive structure since it is assumed
that all measurements are already available. If measurements only occur
sequentially, for example in tracking, then recursive methods such as the
Kalman filter or, for ellipses, the extended Kalman filter are necessary.

## References

  * [1] J. Porrill, “Fitting ellipses and predicting confidence envelopes using a bias corrected kalman filter,” _Image Vision Comput._ , vol. 8, no. 1, pp. 37–41, Feb. 1990. [Online]. Available: http://dx.doi.org/10.1016/0262-8856(90)90054-9
  * [2] I. Frosio and N. Borghese, “Real-time accurate circle fitting with occlusions,” _Pattern Recognition_ , vol. 41, no. 3, pp. 1041 – 1055, 2008, part Special issue: Feature Generation and Machine Learning for Robust Multimodal Biometrics. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0031320307003925
  * [3] S. J. Ahn, W. Rauh, and M. Recknagel, “Ellipse fitting and parameter assessment of circular object targets for robot vision,” in _Intelligent Robots and Systems, 1999. IROS ’99. Proceedings. 1999 IEEE/RSJ International Conference on_ , vol. 1, 1999, pp. 525–530 vol.1.
  * [4] N. Greggio, L. Manfredi, C. Laschi, P. Dario, and M. Carrozza, “Robotcub implementation of real-time least-square fitting of ellipses,” in _Humanoid Robots, 2008. Humanoids 2008. 8th IEEE-RAS International Conference on_ , Dec 2008, pp. 174–181.
  * [5] N. Chernov and C. Lesort, “Statistical efficiency of curve fitting algorithms,” _eprint arXiv:cs/0303015_ , Mar. 2003.
  * [6] Zhengyou Zhang, “Parameter estimation techniques: a tutorial with application to conic fitting,” _Image and Vision Computing_ , vol. 15, no. 1, pp. 59–76, 1997.
  * [7] M. Baum, V. Klumpp, and U. Hanebeck, “A novel bayesian method for fitting a circle to noisy points,” in _Information Fusion (FUSION), 2010 13th Conference on_ , July 2010, pp. 1–6.
  * [8] M. Baum and U. Hanebeck, “Fitting conics to noisy data using stochastic linearization,” in _Intelligent Robots and Systems (IROS), 2011 IEEE/RSJ International Conference on_ , Sept 2011, pp. 2050–2055.