Enhanced English Universal Dependencies: An Improved Representation for Natural Language Understanding Tasks

Many shallow natural language understanding tasks use dependency trees to extract relations between content words. However, strict surface-structure dependency trees tend to follow the linguistic structure of sentences too closely and frequently fail to provide direct relations between content words. To mitigate this problem, the original Stanford Dependencies representation also defines two dependency graph representations which contain additional and augmented relations that explicitly capture otherwise implicit relations between content words. In this paper, we revisit and extend these dependency graph representations in light of the recent Universal Dependencies (UD) initiative and provide a detailed account of an enhanced and an enhanced++ English UD representation. We further present a converter from constituency to basic, i.e., strict surface structure, UD trees, and a converter from basic UD trees to enhanced and enhanced++ English UD graphs. We release both converters as part of Stanford CoreNLP and the Stanford Parser.

Introduction
Since its first version, the Stanford Dependencies (SD) representation (de Marneffe et al., 2006) has had the status of being both a syntactic and a shallow semantic representation. This dual status is also reflected in the usage of SD in natural language processing tasks which broadly fall into two categories. The first category is composed of tasks that require a syntactic tree such as source-side reordering for machine translation (e.g., Genzel (2010)) and sentence compression (Galanis and Androutsopoulos, 2010). For these tasks, a sound syntactic representation is more important than the relations between individual words. The second and much larger category is composed of a wide range of shallow natural language understanding (NLU) tasks such as biomedical text mining (e.g., Airola et al. (2008)), open domain relation extraction (e.g., Mausam et al. (2012)), and unsupervised semantic parsing (Poon and Domingos, 2009). For these tasks, the relations between content words are more important than the overall tree structure. Not surprisingly, we observe a similar divide if we look at which one of the three SD representations is being used for the individual downstream tasks. Most systems that require a syntactic representation use basic SD trees which are guaranteed to be a strict surface syntax tree. On the other hand, most systems that are concerned with the relations between content words use the collapsed or CCprocessed SD representations. These representations may be graphs instead of trees, and may contain additional and augmented relations that explicitly capture otherwise implicit relations between content words. To illustrate one of the differences between these representations, consider the sentence "Fred started to laugh". The basic SD representation of this sentence lacks a direct relation between the controlled verb laugh and its controller, Fred, while in the CCprocessed SD representation, this relation is made explicit with an additional subject edge.
Fred started to laugh nsubj The popularity of these extended representations suggests that their existence plays a major role in the popularity of Stanford Dependencies. In recent years, there has been a lot of interest in extending Stanford Dependencies to other languages, including morphologically rich ones (McDonald et al., 2013;Tsarfaty, 2013;de Marneffe et al., 2014). These individual projects ultimately led to the Universal Dependencies (UD) initiative (Nivre et al., 2016) whose goal is to develop crosslinguistically consistent treebank annotations for as many languages as possible. While this project recognizes the status of dependency formalisms as semantic representations and based many design decisions on their impact on NLU tasks, the majority of efforts so far have focused on the development of the basic UD representation and the annotation of treebanks. Both de Marneffe et al. (2014) and Nivre et al. (2016) also mention an enhanced UD representation and acknowledge its usefulness but neither gives a detailed account of what such a representation should look like. In this paper, we revisit and extend the collapsed and CCprocessed SD representations in light of the recent developments by the Universal Dependencies initiative. We provide a detailed account of an enhanced English UD representation and introduce the enhanced++ representation which we deem even better suited for many NLU tasks. Further, we describe our implementation of a converter from phrase-structure trees to basic UD trees, and a converter from basic to enhanced and enhanced++ English UD graphs. We release both tools as part of Stanford CoreNLP and the Stanford Parser.

The enhanced UD representation
The enhanced English UD representation aims to make implicit relations between content words more explicit by adding relations and augmenting relation names. In the development of this representation, we adhered to the guidelines by Nivre et al. (2016) which state that an enhanced dependency graph may only contain additional dependencies or introduce additional language-specific relations. As a result, enhanced UD graphs contain all the relations of the basic UD tree and the following additional relations. Augmented modifiers One major difference between the SD and UD representations is that the head of prepositional phrases (PP) is the prepositional complement and no longer the preposition itself. Therefore, there already exists a relation between the content word in the prepositional phrase and the word that is being modified by the PP in the basic UD representation, and there is no need for an additional relation in the enhanced representation. However, the collapsed SD graphs of sentences with PPs do not only contain additional relations, they also include the preposition in the relation name. This helps to disambiguate the type of modifier and further facilitates the extraction of relationships between content words, especially if a system incorporates dependency information using very simple methods such as by only considering paths between two nodes. For this reason, all nominal modifiers (nmod) in enhanced UD graphs also include the preposition in their relation name as exemplified in the following phrase. the house on the hill det nmod:on case det The same is true for more complex PPs which are either analyzed as adverbial clause modifiers (advcl) or as adjectival clause modifiers (acl) in UD. If an adverbial clause or an adjectival clause is introduced by a marker (mark), we add the marker to the relation name.
he brushed his teeth after eating dinner nsubj nmod:after mark dobj dobj Augmented conjuncts In a similar manner, enhanced UD graphs also contain conjunct relations that are augmented with their coordinating conjunction. This makes the type of coordination between two phrases more explicit which is particularly useful in phrases with multiple coordinating conjunctions, such as the following phrase.
apples and bananas, or oranges cc conj:and cc conj:or Propagated governors and dependents In basic UD trees of clauses with conjoined phrases, only the first conjunct has explicit relations to the governor and the dependents of the conjoined phrase. In the enhanced UD graph, the implicit relations of the other conjuncts are made explicit with additional relations. In the case of conjoined noun phrases, each noun phrase becomes an argument of the head of the first conjunct, e.g., it becomes the subject of the main verb as in the following example.
Sue and Paul are running Subjects of controlled verbs Basic UD graphs of sentences that contain a controlled embedded verb lack a direct relation between the controlled verb and the controller. Therefore, the enhanced UD graphs contain a special controlling subject (nsubj:xsubj) relation between the embedded verb and the controller.
Sue wants to buy a hat nsubj nsubj:xsubj xcomp mark dobj det

The enhanced++ UD representation
The enhanced representation provides reasonable analyses for most English sentences. However, there are some constructions in English that lead to an analysis which is suboptimal for many NLU systems that try to extract relationships between entities, such as open domain relation extraction (e.g., Mausam et al. (2012)), or extracting relationships between objects in image descriptions (Schuster et al., 2015). One set of problematic constructions involves partitive noun phrases such as both of the girls in which both of the acts semantically as a quantificational determiner. In the basic UD representation, however, both is the head of such a partitive phrase while the semantically very similar phrase both girls is headed by girls: Considering that many relation extraction systems use simple dependency tree patterns to extract entities and their relations, these different analyses are clearly problematic. In the first phrase, the determiner both appears to be the subject while in the second phrase, girls is the subject, but ideally both phrases would be analyzed in a similar way. However, in order to obtain a similar analysis for both phrases, we would have to change the structure of the basic dependency trees, which is not allowed according to the guidelines for enhanced dependency graphs. Another set of problematic constructions involves multiword prepositions such as in front of. As illustrated in the following tree, the enhanced representation of "a house in front of the hill" contains a relation between house and front, and a relation between front and hill.
the house in front of the hill det nmod:in case nmod:of case det But for most tasks, the relation between house and hill is going to be more relevant. This relation could be made explicit by adding another relation between house and hill, but as we are not allowed to delete any relations from the basic representation, we would encode some information twice. These two issues should illustrate that there exist several phenomena for which both the basic UD representation and the enhanced UD representation provide suboptimal analyses. For this reason, we argue for another representation which allows for the deletion of relations from the basic UD tree, as this gives us more flexibility in analyzing several constructions in English, including the ones mentioned above. We therefore introduce enhanced++ UD graphs and propose different analyses for the following common phenomena in English.
Partitives and light noun constructions For the analysis of partitive noun phrases such as both of the girls, we follow Barwise and Cooper (1981) and Keenan and Stavi (1986) and treat the first part of the phrase as a quantificational determiner. We promote the semantically salient noun phrase, e.g., girls in our example, to be the head of the partitive and we analyze the quantificational determiner as a flat multi-word expression that is headed by its first word.
In order to mark that these quantificational determiners are not regular determiners, we attach them using the special relation quantificational modifier (det:qmod).
Both of the girls are reading det:qmod mwe mwe nsubj aux Light noun constructions (Simone and Masini, 2014) such as a panel of experts or a bunch of people pose similar challenges because the light nouns are the head of these phrases in the corresponding basic UD trees. However, just like the partitives, the second noun phrase tends to be the semantically salient one in these constructions while the first part of the phrase again serves as a quantificational determiner. We therefore analyze light noun constructions exactly like partitives as illustrated in the following example. 1 a bunch of people are coming det:qmod mwe mwe nsubj aux Multi-word prepositions As mentioned above, multiword prepositions such as in front of tend to obscure the relation between two content words. While the basic UD representation analyzes some multi-word expressions with function words, e.g., due to, using a special mwe relation, the set of these expressions is very limited and does not include many multi-word prepositions. To introduce a direct relation between content words in the enhanced++ UD representation, we also analyze these multi-word prepositions as flat multi-word expressions headed by the first word, and we attach the head of the phrase to the following noun phrase.
the house in front of the hill det nmod:in front of mwe mwe case det Conjoined prepositions and prepositional phrases Clauses that contain conjoined prepositions such as "I bike to and from work" also pose some challenges. Ideally, the UD graph should encode that there is an nmod:to as well as an nmod:from relation between bike and work. Further, we also want to encode that bike to work and bike from work are conjoined by and. In order to be able to represent all of this information, the CCprocessed SD representation introduced copy nodes which we adapt in the enhanced++ representation. The analysis of this example then contains a copy of bike, namely bike . This copy node is attached to the original node as a conjunct resulting in the following UD graph.
I bike bike to and from work nmod:to nmod:from case cc conj:and nsubj nsubj conj:and Note that this graph contains the same relations between content words as the basic UD tree of the clause "I bike to work and I bike from work". Similar complexities arise with clauses that contain a conjoined prepositional phrase such as "She flew to Bali or to Turkey". Again, the UD graph should encode that there is an nmod:to relation between flew and Bali and another nmod:to relation between flew and Turkey and at the same time it should encode that these two relations are conjoined by or. For these reasons, we also analyze such clauses with copy nodes. Relative pronouns We also analyze relative pronouns differently in the enhanced++ representation as compared to the basic UD representation. Similar to the collapsed SD representation, the referent of the pronoun is directly attached to the governor of the pronoun. Further, we attach the relative pronoun to its referent with a referent (ref) relation. The following example illustrates the differences between basic UD and enhanced++ UD graphs which contain relative clauses. In the basic UD tree (a), the head of the relative pronoun who is lived and there is no direct relation indicating that boy is an argument of lived. In the enhanced++ representation (b), on the other hand, boy is the subject of lived, and who is attached to its referent, boy.

Generating dependency trees and graphs
There exist multiple ways to obtain trees and graphs in the various UD representations for a given sentence. The basic UD trees can be either generated directly by using a dependency parser or by using a constituency parser followed by a converter from phrase structure to dependency trees. The enhanced and enhanced++ dependency graphs can be obtained by post-processing basic dependency trees. In the following two sections, we describe our converter from phrase structure to basic English UD trees and how to obtain enhanced and enhanced++ dependency graphs from basic English UD trees.

Converting to basic dependencies
Our converter from phrase structure to dependency trees is based on the Stanford Dependencies converter (de Marneffe et al., 2006) which we updated according to the English Universal Dependencies guidelines.
To determine the structure of the dependency tree, we use a semantic head finder which operates similarly to the Collins head finder (Collins, 1999). For each constituent type, we define a set of rules that determine from which of its children the constituent inherits its head; terminal nodes have their surface token as head. These rules are mostly conditions on the constituent type of the child but unlike the classic Collins head finder rules, some of them also take surface tokens into account which is necessary for distinguishing between main verbs and auxiliaries. We traverse the constituency tree in depth-first order and use these rules to obtain and store the head of each constituent, resulting in a tree in which every node has exactly one surface token as head. The head of each token is then simply the head of its lowest ancestor whose head is not the token itself.
To determine the relation types, we define for each grammatical relation a set of tree patterns in the form of tregex expressions (Levy and Galen, 2006). For each headdependent pair we try to find a pattern that matches the subtree rooted at the lowest common ancestor of the head and the dependent. If such a pattern exists, we assign its corresponding grammatical relation to the head-dependent pair.
In the rare cases where no pattern matches, we assign the most general relation dep. This procedure allows us to obtain correct dependency trees in most cases. Two phenomena, however, require additional consideration. First, as previously mentioned, the UD representation defines several multi-word expressions with function words that behave like a single word such as because of or in case. Extracting the correct structure for these expressions is often challenging because many of these expressions are not a constituent according to the Penn Treebank annotation guidelines (Marcus et al., 1993). We resolve this issue by preprocessing all phrase structure trees that contain multi-word expressions such that the entire expression forms a constituent. Second, the outlined procedure often attaches wh-words in questions to the wrong head. For a question such as What does Peter seem to have?, our procedure would attach what to the head of the matrix clause, seem, instead of the head of the embedded clause, have. If we were only concerned with converting manually annotated treebanks, we could resolve these ambiguities by making use of the indexed empty nodes in the phrase structure trees, as proposed by Choi and Palmer (2012). However, the output of most constituency parsers does not contain these empty nodes. Therefore, we try to solve this issue by considering the selectional restric-tions of the verb in the matrix clause and if the attachment of the wh-word violates these restrictions we try to reattach it to the head of the embedded clause.

Converting to enhanced and enhanced++ dependencies
For the conversion from basic dependencies to enhanced and enhanced++ dependencies we mainly rely on dependency tree patterns in the form of Semgrex expressions (Chambers et al., 2007). Semgrex expressions allow one to match subgraphs of a dependency graph based on properties of the nodes and their relations. For most of the enhancements, we search for syntactic patterns and then either modify the relation name or introduce new relations. The multi-word prepositions and the quantificational constructions, however, do not match any distinct dependency patterns. For these modifications, we rely on lists of specific expressions and modify the graph structure whenever we encounter one of them. Our converter generates enhanced and enhanced++ dependency graphs as described in the previous sections, with one exception. Currently, we don't propagate object or nominal modifier relations in clauses with conjoined verb phrases such as "the store buys and sells cameras". The reason for this is that there are also many cases such as "she was reading or watching a movie" where it would be wrong to add these relations and there are no syntactic cues that would allow us to distinguish these cases. Nyblom et al. (2013) successfully used a machine learning approach to solve this problem for Finnish but as there currently exists no corpus annotated with enhanced English UD graphs, we leave this to future work.

Evaluation
We evaluate our basic converter against the manually checked English UD treebank v1.1 (Nivre et al., 2015) which contains annotations for all sentences in the EWT corpus (English Web Treebank, Linguistic Data Consortium release LDC2012T13). We convert all phrase structure trees in the EWT corpus to basic UD trees and compare the output with the manually checked trees using the official CoNLL Shared Task evaluation script. The results of this evaluation are presented in Table 1.
These results indicate that our converter is able to convert phrase structure trees from a variety of genres to basic UD trees at high accuracy. If we compare the performance of the converter across the individual genres, we can see that the converter performs best on sentences from weblogs and newsgroups and slightly worse on sentences from emails, from an online question-answering site, and from online business reviews. A qualitative error analysis showed that the main reason for the small drop in performance on the question-answer and review corpora is that these corpora contain a lot of ungrammatical sentences written by nonnative English speakers. The main reason for the lower performance of the converter on the email corpus is that this corpus contains a lot of corporate email signatures whose corresponding phrase structure trees consist of a single flat fragment from which it is very hard to extract properly structured dependency trees.  Table 1: Labeled attachment score (LAS), unlabeled attachment score (UAS), and accuracy of the converter from phrase structure trees to basic English UD trees on the individual genres of the English Web Treebank corpus. We use the manually corrected English UD corpus v1.1 (Nivre et al., 2015) as a gold standard.

Applications
We believe there are two main applications of our converters. First, our basic UD converter can be used to automatically convert existing treebanks of phrase structure trees to treebanks of UD trees for training dependency parsers. For example, we have successfully converted the entire Penn Treebank (Marcus et al., 1993) to train models for a neural network dependency parser (Chen and Manning, 2014). Second, our converters can be used either in combination with a constituency parser or a dependency parser to obtain UD graphs for any sentence which can then be utilized in downstream NLU tasks. Schuster et al. (2015), for example, used a preliminary version of the converter to obtain enhanced++ UD graphs from constituency trees. This system uses the UD graphs as input for a parser from image descriptions to a scene representation that captures relationships between objects in a visual scene. Further, the open domain information extractor in CoreNLP (Angeli et al., 2015) also already uses UD graphs to extract relations between entities.

Comparison to AMR
Representing the meaning of sentences as directed graphs has a long tradition in computational linguistics, which goes back to at least Shieber (1984). One graph-based semantic representation that received significant attention in recent years is the Abstract Meaning Representation (AMR) (Banarescu et al., 2013). AMR also encodes sentences as directed graphs but compared to UD graphs, it aims to abstract further away from the surface form of sentences. To achieve this goal, it encodes sentences using PropBank framesets (Palmer et al., 2005) and approximately 100 fixed relations. This makes AMR a deeper and more canonicalized semantic representation as compared to UD graphs. While these are obviously desirable properties, we nevertheless believe that our representation has some advantages over AMR, especially when it is being used in shallow natural language understanding tasks. In terms of expressivity, UD graphs have the advantage that they encode the meaning of sentences in terms of relations between surface form tokens and they are therefore as expressive as natural language. The expressivity of AMR, on the other hand, is constrained by the finite set of PropBank framesets. This is particularly problematic for neologisms such as to venmo for which no corresponding PropBank framesets exist. Further, we are using existing resources of languages to disambiguate argument types, e.g., by including the preposition in the relation name of nominal modifiers, which avoids many hard labeling decisions for parsers as well as for human annotators. While this approach is occasionally too simplistic as our representation does not, for example, distinguish between comitative with and instrumental with, it nevertheless sufficiently disambiguates argument types for most domain-specific NLU tasks. AMR, on the other hand, requires the labeling of every argument with a semantic role which is -apart from labeling agents and patientsa very hard task. Another appeal of UD graphs is their simplicity. A lot of the frequent relations, such as nominal subject or object are very intuitive which makes them easily interpretable by non-experts. Compared to UD graphs, AMR graphs are more complex and require more explanation even if one is only interested in extracting simple relations such as subject-verb-object triplets. Finally, from a practical point of view, sentences can be parsed to UD graphs with much higher accuracy than to AMR graphs. While the existence of high-performing parsers clearly should not be the main desideratum in the design of a semantic representation, this aspect plays ultimately an important role in the usefulness of a representation.

Limitations
As explained above, most of the additional relations in the enhanced and enhanced++ representations can be added with syntactic rules. This purely syntactic approach tends to work well in practice, and has the appeal of being very simple and easily comprehensible. However, in some cases this approach leads to UD graphs that encode a different meaning than the original sentence, which can be problematic in downstream NLU tasks. The issue with this UD graph is that it does not encode the meaning of the original sentence but instead encodes approximately the meaning of the sentence "Everybody wants that everybody buys a house". In order to preserve the original meaning and to encode the relation between the controlled verb and its subject we would have to introduce variables and consequently abandon one of our core principles, namely that we encode the meaning of a sentence in terms of relations between surface form tokens.
2 Thanks to Christopher Potts for pointing out this issue.
A second issue concerns the propagation of dependents. Note that in the case of conjoined verbs or verb phrases, we are effectively performing a reverse conjunction reduction which can lead to problematic analyses in combination with generalized quantifiers. For example, consider the sentence "Everybody sleeps or is awake" with the following enhanced UD graph.
Everybody sleeps or is awake nsubj nsubj conj:or cop cc The issue with this UD graph is that it approximately encodes the meaning of the sentence "Everybody sleeps or everybody is awake" which again differs from the meaning of the original sentence. Lastly, another issue concerns sentences with conjoined subjects, such as "Sue and Mary are carrying a piano".
Sue and Mary are carrying a piano nsubj nsubj conj:and cc dobj det aux Unlike in the previous two examples, the issue with this UD graph is not that it encodes a different meaning than the original sentence, but rather that this representation of conjoined subjects favors a distributive interpretation. Ideally, we would propagate the subject relation only when it is clear that a distributive interpretation is intended, but as the question whether a conjoined subject should be interpreted distributively or collectively also depends on world knowledge and the context, we are not able to make this distinction based on the information that is contained in a single sentence.

Conclusion
In this paper, we presented the first detailed account of an enhanced English Universal Dependencies representation. We further argued for additional modifications of the tree structure to facilitate extracting relations between content words and described these modifications as part of the enhanced++ representation. Finally, we described how both of these representations can be automatically generated from phrase structure trees or basic dependency trees with high accuracy.

Acknowledgements
This work greatly benefitted from many improvements to the original Stanford Dependencies converter by Marie-Catherine de Marneffe and John Bauer, which we gratefully acknowledge. We also thank Natalia Silveira and Timothy Dozat for helpful discussions, and the entire Universal Dependencies community for developing the guidelines and their valuable input. This work was supported in part by a gift from IPSoft, Inc. The first author is also supported by a Goodan Family Graduate Fellowship.


extfieldsofstudy:
Mathematics
Computer Science

text:
Graph clustering, variational image segmentation methods and Hough transform scale detection for object measurement in images

We consider the problem of scale detection in images where a region of interest is present together with a measurement tool (e.g. a ruler). For the segmentation part, we focus on the graph based method by Flenner and Bertozzi which reinterprets classical continuous Ginzburg-Landau minimisation models in a totally discrete framework. To overcome the numerical difficulties due to the large size of the images considered we use matrix completion and splitting techniques. The scale on the measurement tool is detected via a Hough transform based algorithm. The method is then applied to some measurement tasks arising in real-world applications such as zoology, medicine and archaeology.

Introduction
Image segmentation denotes the task of partitioning an image in its constituent parts. Feature based segmentation looks at distinctive characteristics (features) in the image, grouping similar pixels into clusters which are meaningful for the application at hand. Typical examples of features are based on greyscale/RGB intensity and texture. Mathematical methods for image segmentation are mainly formalised in terms of variational problems in which the segmented image is a minimiser of an energy. The most common image feature encoded in such energies is the magnitude of the image gradient, detecting regions (or contours) where sharp variations of the intensity values occur. Examples include the Mumford-Shah segmentation approach [49], the snakes and geodesic active contour models [39,19]. Moreover, in [21] Chan and Vese proposed an instance of the Mumford-Shah model for piecewise constant images whose energy is based on the mean greyvalues of the image inside and outside of the segmented region rather than the image gradient and hence does not require strong edges for segmentation. Other image segmentation methods have been considered in [40,29]. They rely on the use of the total variation (TV) seminorm [5], which is commonly used for image processing tasks due to its properties of simultaneous edge preservation and smoothing (see [55]).
The non-smoothness of most of the segmentation energies renders their numerical minimisation usually difficult. In the case of the Mumford-Shah segmentation model the numerical realisation is additionally complicated by its dependency on the image function as well as the object contour. To overcome this, several regularisation methods and approximations have been proposed in the literature, e.g. [4,12,13,66] for Mumford-Shah segmentation. In the context of TV-based segmentation models the Ginzburg-Landau functional has an important role. Originally considered for the modelling of physical phenomena such as phase transition and phase separation (cf. [14] for a survey on the topics) it is used in imaging for approximating the TV energy. Some examples of the use of this functional in the context of image processing are [29,27,28], which relate to previous works by Ambrosio and Tortorelli on diffuse interface approximation models [5,4].
Such variational methods for image segmentation have been extensively studied from an analytical point of view and the segmentation is usually robust and computationally efficient. However, variational image segmentation as described above still faces many problems in the presence of low contrast and the absence of clear boundaries separating regions. Their main drawback is that they are limited to image features which can be mathematically formalised (e.g. in terms of an image gradient) and encoded within a segmentation energy. In recent years dictionary based methods have become more and more popular in the image processing community, complementing more classical variational segmentation methods. By learning the distinctive features of the region to be segmented from examples provided by the user, these methods are able to segment the desired regions in the image correctly.
In this work, we consider the method proposed in [9,10,11,43] for image segmentation and labelling. This approach goes beyond the standard variational approach in two respects. Firstly, the model is set up in the purely discrete framework of graphs. This is rather unusual for variational models where one normally considers functionals and function spaces defined on subdomains of R 2 in order to exploit properties and tools from convex and functional analysis and calculus of variations. Secondly, the new framework allows for more flexibility in terms of the features considered. Additional features like texture, light intensity or others, can be considered as well without encoding them in the function space or the regularity of the functions. Due to the possibly very large size of the image (nowadays of the order of megapixel for professional cameras) and the large number of features considered, the construction of the problem may be computationally expensive and often requires reduction techniques [51,50,30]. In several papers (see, e.g., [59,61,35]) this problem was approached by means of the graph cut objective function.
In this paper, we also address the problem of detection of objects with geometrical properties that are a priori known. An example is the detection of lines and circles. These objects can be identified by mapping them onto an auxiliary space where relevant geometrical properties (such as linear alignment and roundness) are represented as peaks of specific auxiliary functions. In this work, we use the Hough transform [37] to detect measurement tools (rulers, concentric circles of fixed radii) with the intent of providing quantitative, scale-independent measurements of the region segmented by one of the techniques described above. In this way, an absolute measurement of the region of interest in the image is possible, independent of the scale of the image, which could depend, for instance, on the distance of the objective to the camera.
We demonstrate the use of our method in the context of real world applications in which segmentation and subsequent object measurement are crucial. Our main application is the measurement of the white forehead patch (blaze) of male pied flycatchers, which has been studied with regard to sexual selection in [58], see Figure 1. The forehead patch is known to vary between individuals [41] and can be subject to both intra- [38] and intersexual [52] selection with pied flycatchers from Spain preferring males with large patches. Forehead patch size has been shown to signal male phenotypic quality through plasma oxidative stress and antioxidant capacity [48]. However, in all studies to date the measurements of patches have been inconsistent and generally inaccurate. For example some studies have simply measured patch height [25], whereas Potti and Montalvo [52] assumed the shape to be a trapezium with area equal to 0.5(B + b)H, B being the white patch width, b the bill width and H the height of the white patch. Morales et al. [46] measured the length and breadth of the forehead patch with callipers to the nearest 0.01mm and its size (mm 2 ) was calculated as the area of a rectangle. Other studies have measured the patches from photographs, e.g., Järvistö et al. [38] Ruuskanen et al. [56] and Sirkiä et al. [60] who photographed the forehead with a scale bar included in each picture, and measured the patch as the white area in mm 2 using IMAGEJ software [2]. But none of these three papers provide methods of how the measurement was actually achieved, e.g., whether patches were delineated or roughly estimated with a simple shape. Most recently Moreno et al. [48] analysed digital photos of forehead patches with Adobe PhotoShop CS version 11.0. relating the distance of 1mm on the ruler to number of pixels, and used this to estimate length. Zooming to 400% and using the paintbrush tool with 100% hardness and 25% spacing the authors delineate the patch and measure the area of the white areas on forehead. While this is the best measurement method to date, it still is subject to human measurement error and subjective assessment of patch boundaries. We report some segmentation results obtained by manual selection and polygon fitting in Figure 2. In this manuscript we use a mathematically robust approach to segment the blaze independently to provide an accurate measurement of forehead patch area. A similar challenge can be encountered in medical applications monitoring and quantifying the evolution of skin moles for early diagnosis of melanoma (skin cancer). A normally user-dependent measurement of the mole is performed using a ruler located next to it. A picture is then taken and used for future comparisons and follow-  Figure 2: Flycatcher blaze segmentation of the images 1b and 1c obtained either by using the 'magic wand' tool of the IMAGEJ software, similarly as described by Moreno [48] or by trapezium fitting as suggested by Potti and Montalvo [52]. In the first case the result is strongly user-dependent, in the second one the blaze area is overestimated.
up, see Figure 3 and compare [20,1] for previous attempts of automatic detection of melanomas. For such an application, a systematic quantitative analysis is also required * .
In several other applications the task of measuring objects directly from the image is encountered. These include zoological and behavioural studies arising in the animal world where detecting size, shape and possible symmetries of specific distinctive animal features can be useful, as well as, for instance, in archaeological digs where the measurement of finds is important for comparisons and classification [36]. Figure 3: The monitoring and measuring of moles is essential for the early diagnosis of melanoma. Normally, due to their small size, they can be measured by juxtaposing a small ruler with them.
Outline of the method. We consider the image as a graph whose vertices are the image pixels. Similarity between pixels in terms of colour or texture features is modelled by a weight function defined on the set of vertices. Our method runs as follows. Firstly, using examples provided by the user (dictionaries) as well as matrix completion and operator splitting techniques, the segmentation of the region of interest is performed. In the graph framework, this corresponds to cluster together pixels having similar features. This is obtained by minimising on the graph the Ginzburg-Landau functional typically used in the continuum setting to describe diffuse interface problems. In order to provide quantitative measurements of the segmented region, a second detection step is then performed. The detection here aims to identify the distinctive geometrical features of the measurement tool (such as line alignment for rulers or circularity for circles) to get the scale on the measurement tool considered. The segmented region of interest can now be measured by simple comparisons and quantitative measurements such as perimeter and area can be provided.

Contribution
We propose a self-contained programme combining automated detection and subsequent size measurement of objects in images where a measurement tool is present. Our approach is based on two powerful image analysis techniques in the literature: a graph segmentation approach which uses a discretised Ginzburg-Landau energy [9] for the detection of the object of interest and the Hough transform [37] for detecting the scale of the measurement tool. While these methods are state of the art, their combination for measuring object size in images proposed in this paper is new. Moreover, to our knowledge there is only little contribution in the literature that broach the issue of how the graph segmentation approach as well as the Hough transform are applied to specific problems [10,11,34]. Indeed, here we present these methodologies in detail, especially discussing important aspects in their practical implementation, and demonstrate the robust applicability of our programme for measuring the size of objects, showcasing its performance on several examples arising in zoology, medicine and archaeology. Namely, we first apply our combined model for the measurement of the blaze on the forehead of male pied flycatchers, for which we run a statistical analysis on the accuracy and predicted error in the measurement on a database of thirty images. State-of-the-art methods for such a task typically require the user to fit polygons inside or outside the blaze [52] or to segment the blaze by hand [48]. Similarly, the scale on the measurement tool is typically read from the image by manually measuring it on the ruler. With respect to medical applications, we apply our combined method for the segmentation and measurement of melanomas. Although efficient segmentation methods for automatised melanoma detection already exist in literature (see, e.g., [20,1]), up to the knowledge of the authors no previous methods providing their measurement by detecting the scale on the the ruler placed next to them (see Figure 3) exist. Conversely, in the case of archaeological applications, some models for the automatic detection of the measurement tool in the image exist [36] but no automatic methods are proposed for the segmentation of the region of interests.
Organisation of the paper. In Section 2 we present the mathematical ingredients used for the design of the graph based segmentation technique used in [9,10,11,43]. They come from two different worlds: the framework of diffusion PDEs used for modelling phase transition/separation problems (see Section 2.1) and graph theory and clustering, see Section 2.2. In view of a detailed numerical explanation, we also recall a splitting technique and a popular matrix completion technique used in our problem to overcome the computational costs. In Section 3 we review the Hough transform and explain how that is applied to the scale detection problem. Finally, Section 4 contains the numerical results obtained with our combined method applied to the problems described above. For completeness, we report in Appendix A some variational models for image segmentation mentioned in the Introduction.

Image segmentation as graph clustering
We present in this section the mathematical background for the design of the Ginzburg-Landau (GL) graph segmentation algorithm introduced in [9]. There, the image segmentation problem is rephrased as a minimisation problem on a graph defined by features computed from the image. Compared to the methods above, the graph framework allows for more freedom in terms of the possible features used to describe the image, such as texture.

The Ginzburg-Landau functional as approximation of TV
In the following, we recall the main properties of the original continuum version of the GL functional explaining its importance in the context of image segmentation problems as well as the main concepts of graph theory which will be used for the segmentation modelling.
Several physical problems modelling phase transition and phase separation phenomena are built around the well-known GL functional: The functional above is defined in the continuous setting. Here, Ω represents a open subset of R d , d = 2, 3, u : Ω → R is the density of a two-phase material and W (u) is a double-well potential, e.g. W (u) = 1 4 (u 2 − 1) 2 . The two wells ±1 of W correspond to the two phases of the material. The parameter ε > 0 is the spatial scale. Variational models built around this functional are also referred to as diffuse interface models because of the interface appearing between the two regions containing the phases (i.e. the two wells of W ) due to the competition between the two terms of the functional (2.1). Nonetheless, some smoothness preventing u from having jumps between the two regions is ensured by the first regularisation term.
The use of the GL functional has become very popular in image processing due to its connections with the total variation (TV) seminorm. In [44,45], for instance, Γ-convergence properties of (2.1) to the TV functional are shown. Thus, the GL functional is very often used as a quadratic approximation of total variation. Fast numerical schemes relying on these connections have been designed for many imaging problems, thus overcoming the issues related to nonsmooth TV minimisation [5,29,21]. In image processing, the functional considered often is of the form E(u) := GL(u) + λ φ(u, u 0 ), (2.2) where φ(u, u 0 ) is a fidelity term measuring the distance of the reconstructed image u to the given image u 0 . Depending on the application, different data fidelities are employed. Typically, they are related to statistical and physical assumptions of the model considered. Standard examples of fidelity terms are φ(u, u 0 ) = u − u 0 d L d (Ω) , d = 1, 2. The parameter λ > 0 determines the influence of the data fit compared to the regularisation. Taking the L 2 gradient descent of (2.2) we get the following evolutionary PDE, known in the literature as the Allen-Cahn equation [3] with an additional forcing term due to the fidelity φ: Steady states of equation (2.3) are the stationary points of the energy E in (2.2). Note that E is not convex so uniqueness is not guaranteed and, consequently, the long time behaviour for solutions of (2.3) will depend on the initial condition. The linear diffusion term weighted by ε appearing in (2.3) allows for fast solvers using for instance the Fast Fourier Transform (FFT) which translates the Laplace operator into a multiplication operator on the Fourier modes.

Towards the modelling: the graph framework
In the following, we rely on the method presented in [9,43] for high-dimensional data classification on graphs which has been applied to several imaging problems [10,11], showing good performance and robustness. We consider the problem of binary image segmentation where we want to partition a given image into two components where each component is a set of pixels (also called a cluster, or a class) and represents a certain object or group of objects. Typically, some a priori information describing the object(s) we want to extract is given and serves as initial input for the segmentation algorithm. For image labelling, in [9] two images are taken as input: the first one has been manually segmented in two classes and the objective is to automatically segment the second image using the information provided by the segmentation of the first one. We revise in the following the main ingredients of the model considered and start from a quick review of concepts in graph theory. We consider an image I as a collection of pixels arranged in a rectangle of size S := N × M , i.e.
For each x ∈ I, we define the image neighbourhood of x as the set N(x) := {y ∈ I : |x 1 − y 1 | ≤ τ and |x 2 − y 2 | ≤ τ }, with τ ∈ N fixed, i.e. N(x) contains the pixels in a (2τ + 1) × (2τ + 1) sized square centred at x. For some appropriate K ∈ N, we associate to every vertex x ∈ I a vector z ∈ R K encoding selected characteristics of the neighbourhood N(x). These characteristics are related to the grey or RGB (red, green, blue) intensity values as well as the texture features of the neighbourhood. In Section 2.5, we will explain in more detail our feature vector construction. The map ψ : I → R K , x → z is called the feature function. With an eye to constructing the feature vectors, it will be useful to associate a neighbourhood vector ν(x) := (x j ) j∈N(x) ∈ I (2τ +1)×(2τ +1) to each neighbourhood, such that the ordering of the x j in ν(x) is consistent between pixels x, e.g., order the pixels from each square N(x) from left to right and top to bottom. The specific choice of ordering is not important, as long as it is consistent for each pixel neighbourhood.
Next we construct a simple weighted undirected graph G = (V, E, w) whose vertices correspond to the pixels in I and with edges whose weights depend on the feature function ψ. Let V be a vertex set of cardinality S. To emphasize that each vertex in V corresponds to exactly one pixel in I, we will label the vertex corresponding to x ∈ I by x as well. Let w : V × V → R be a symmetric and nonnegative function, i.e.
We define the edge set E as the collection of all undirected edges connecting nodes x i and x j for which w(x i , w j ) > 0 [23]. The function w restricted to E ⊂ V × V is then a positive edge weight function. In our applications we define w as whereŵ : R K × R K → R is a given function and ψ is the feature function.
In operator form, the weight matrix W ∈ R S×S is the nonnegative symmetric matrix whose elements are w i,j = w(x i , x j ). In the following, we will not distinguish between the two functions w andŵ and, with a little abuse of notation, we will write w(z i , z j ) forŵ(z i , z j ).
Remark 2.1. Weight functions express the similarities between vertices and will be used in the following to partition V into clusters such that the sum of the edge weights between the clusters is small. There are many different mathematical approaches to attempt this partitioning. When formulated as a balanced cut minimisation, the problem is NP-complete [67], which inspired relaxations which are more amenable to computational approaches, many of which are closely related to spectral graph theory [59]. We refer the reader to [23] for a monograph on the topic. The method we use in this paper can be understood (at least in spirit, if not technically, [63,64]) as a nonlinear extension of the linear relaxed problems.
To solve the segmentation problem, we minimise a discrete GL functional (which is formulated in the graph setting, instead of the continuum setting), via a gradient descent method similar to the one described in Section 2.1. In particular, in this setting the Laplacian in (2.3) will be a (negative) normalised graph Laplacian. We will use the spectral decomposition of u with respect to the eigenfunctions of this Laplacian. In Section 2.4 we discuss the Nyström method, which allows us to quickly compute this decomposition, but first we introduce the graph Laplacian and graph GL functional.
The discrete operators. We start from the definition of the differential operators in the graph framework.
For each vertex x ∈ V , we define the degree of x, In operator form, the diagonal degree matrix D ∈ R S×S is defined to have diagonal elements d i,i = d(x i ). A subset A of the vertex set V is connected if any two vertices in A can be connected by a path (i.e. a sequence of vertices such that subsequent vertices are connected by an edge in E) such that all the vertices of the path are in A. A finite family of sets A 1 , . . . , A t is called a partition of the graph if A i ∩ A j = ∅ for i = j and We now have all the ingredients to define the graph Laplacian. Denoting by V the space of all the functions V → R, the graph Laplacian is the operator L : V → V such that: (2.5) We are considering a finite graph of size S, so real valued functions can be identified as vectors in R S . We can then write the graph Laplacian in matrix form as L = D − W or element-wise as: It is worth mentioning (see Remark 2.2 below) that this graph Laplacian is a positive semidefinite operator. Note that by convention the sign of the discrete Laplacian is opposite to that of the (negative semidefinite) continuum Laplacian. The associated quadratic form of L is The quadratic form Q can be interpreted as the energy whose optimality condition corresponds to the vanishing of the graph Laplacian in (2.6).
The eigenvector corresponding to λ 1 is the constant S-dimensional vector 1 S , see [67].
The operator in (2.5)-(2.6) is not the only graph Laplacian appearing in the literature. To set it apart from others, it is also referred to as the unnormalised or combinatorial graph Laplacian. Such operator can be related to the standard continuous differential one through nonlocal calculus [32]. More precisely, the eigenvectors of L converge to the eigenvectors of the standard Laplacian, but in the large sample size limit a proper scaling of L is needed in order to guarantee stability of convergence to the continuum operator [9,11]. Hence, we consider in the following the normalisation of L given by the symmetric graph Laplacian Clearly, the matrix L s is symmetric. Other normalisations of L are possible, such as the random walk graph Laplacian (see [23,67,64]).
In [59, Section 5] a quick review on the connections between the use of the symmetric graph Laplacian (2.8) and spectral graph theory is given. Computing the eigenvalues of the normalised symmetric Laplacian corresponds to the computation of the generalised eigenvalues used to compute normalised graph cuts in a way that the standard graph Laplacian may fail to do, compare [23]. Typically, spectral clustering algorithms for binary segmentation base the partition of a connected graph on the eigenvector corresponding to the second eigenvalue of the normalised Laplacian, using, for example, k-means. For further details and a comparison with other methods we refer the reader to [59] and to [9, Section 2.3] where a detailed explanation on the importance of the normalisation of the Laplacian is given.
The discrete GL functional. Recalling (2.1)-(2.2) and (2.7), we define the discrete GL functional † as Here u 0 represents known training data provided by the user. As before, W (u(x)) = 1 4 (u 2 (x) − 1) 2 is the doublewell potential. The function χ : V → {0, 1} is the characteristic function of the subset of labelled vertices V lab ⊂ V , i.e. χ = 1 on V lab and χ = 0 on V unlab := V c lab . Hence, the corresponding fidelity term enforces the fitting between u and u 0 in correspondence of the the known labels on the set V lab , while the labelling for the pixels in V unlab is driven by the first two regularising terms in (2.9). † 'Discrete GL functional with a data fidelity term' would be a more accurate name, but we opt for brevity here.
The corresponding 2 gradient flow for (2.9) reads The idea is to design a semi-supervised learning (SSL) approach where a priori information for the set V lab (i.e. cluster labels) is used to label the points in the set V unlab . The comparison uses the weight function defined in (2.4) to build the graph by comparing the feature vectors at each point.
Remark 2.3 (The choice of the weight function). As pointed out in [9, Section 2.5], the main criteria driving the choice of the weight function are the desired outcome and the computational efforts required to diagonalise the corresponding matrix W . A common weight function is the Gaussian function Note that this function is symmetric: w(x, y) = w(y, x).
Several approaches to SSL using graph theory have been considered in literature, compare [24,32]. The approach presented here adapts fast algorithms available for the efficient minimisation of the continuous GL functional to the minimisation of the discrete one in (2.9) . In particular, to overcome the high computational costs, we present in the following an operator splitting scheme and a matrix completion technique applied to our problem.

Convex splitting
Splitting methods are used in the study of PDEs. Here, we focus on convex splitting, which is used to numerically solve problems with a general gradient flow structure of the form Here, Ω ⊂ R d is a regular domain, H is the energy functional to minimise and ∇ V indicates formally the Fréchet derivative with respect to the metric in a Banach space V . The main idea of convex splitting is to split the functional H into a convex and a concave part: where both H 1 and H 2 are convex for all u. Denoting by U n the numerical approximation of u(·, n∆t), ∆t > 0, n ≥ 0, a semi-implicit discretization of (2.11)-(2.12) reads The advantage of this method consists in treating the convex part implicitly in time and the concave part explicitly. Typically, nonlinearities are considered in the explicit part of the splitting and their instability is balanced by the effect of the implicit terms. Following [9, Section 3.1], the splitting choice for the discrete GL functional (2.9) is based on the convex splitting of its continuous version with fidelity term, cf. [8,57], and reads where the constant C > 0 has to be chosen large enough such that GL d,2 is convex for u around the wells of W . The differential operator contained in the implicit component of the splitting, GL d,1 , is the symmetric graph Laplacian, which can be diagonalised quickly and inverted using Fourier transform methods. In [9, Section 3.1], more details of the splitting are presented. The time-discretised scheme (2.13) now reads Towards the numerical realisation. The numerical strategy we intend to use is based on the following steps (see Section 2.5 for more details): • At each time step n∆t, n ≥ 0, consider at every point the spectral decomposition of U n with respect to the eigenvectors v k of the operator L s as with coefficients α n . Similarly, use spectral decomposition in the {v k } basis for the other nonlinear quantities appearing in (2.15).
• Having fixed the basis of eigenfunctions, the numerical approximation in the next time step U n+1 is computed by determining the new coefficients α k n+1 in (2.16) for every k through convex splitting (2.15).
The only possible bottleneck of this strategy is the computation of the eigenvectors v k of the operator L s , which, in practice, can be computationally costly for large and non-sparse matrices W . To mitigate this potential problem, we use the Nyström extension (Section 2.4).

Matrix completion via Nyström extension
Following the detailed discussion in [9, Section 3.2], we present here the Nyström technique for matrix completion [51] used in previous works by the graph theory community [30,6] and applied later to several imaging problems [50]. In our problem, the Nyström extension is used to find an approximation of the eigenvectors v k of the operator L s . We will freely switch between the representation of eigenvectors (or eigenfunctions) as a real-valued functions on the vertex set V and as a vectors in R S .
Consider a fully connected graph with vertices V and the set of corresponding feature vectors ψ Thus, finding the spectral decomposition of L s boils down to diagonalising the operator D −1/2 W D −1/2 . This is not obviously easier, as the matrix W , despite being nonnegative and symmetric, may be large and non-sparse, so the computation of its spectrum may be computationally hard. Here, however, we take advantage of the Nyström extension. Given the eigenvalue problem for every point y ∈ V , we approximate the sum on the left hand side using a standard quadrature rule where the interpolation points are chosen by randomly selecting a subset of L points from the set V and the interpolation weights are chosen correspondingly. The Nyström extension for (2.18) then where With this method we can approximate the values of an eigenvector v of W , corresponding to the eigenvalue θ, in the whole set of points V using its values in the subset X and solving the interpolated eigenvalue equation above. Generally, this is not as immediate as it sounds since the eigenvectors of W are not known in advance, however, by choosing y = x j , j = 1, . . . , L, in (2.20), we find an eigenvalue problem for the known matrix with entries w(x j , x i ), which is a much smaller matrix than the full matrix W : If L is small enough such that this eigenvalue problem can be solved, then θ and v(x i ), i = 1, . . . , L, can be computed, which in turn can be substituted back into (2.20) to find an approximation to v(y), for any y ∈ V . In short, we approximate the eigenvectors in (2.18) by extensions of the eigenvectors in (2.21), using the extension equation (2.20), and we approximate the eigenvalues in (2.18) by the eigenvalues from (2.21). The main Nyström assumption is that these approximated eigenvectors and eigenvalues approximately diagonalise W . For numerical implementation and further clarification of the method, it is useful to write it in matrix language. Let us define first the sub matrices W XX ∈ R L × R L and W XY ∈ R L × R S−L as Analogous definitions hold for W Y Y and W Y X . Each of these matrices represents the sub matrix having as elements the weights between the points in X, Y or between the two sets. With this notation, the whole matrix W ∈ R S ×R S can be written in block-form as Similarly, vectors v ∈ R S can be written as We focus on the spectral decomposition of the first block of W , that is W XX . Since this matrix is symmetric, calling Θ X the matrix Θ X = diag(θ 1 , . . . , θ L ) containining the eigenvalues of W XX , then by the spectral theorem W XX = V X Θ X V T X (compare with (2.21)), with V X be the orthogonal matrix having as columns the eigenvectors of W XX . Writing (2.20) for y ∈ Y , in operator form, we The approximated eigenvectors of W can then be written in matrix form as (2.23) Let us observe that Therefore, Nyström extension can be interpreted as the approximation The quality of the approximation of the full W is quantified by the norm of the Schur complement W Y Y − W Y X W −1 XX W XY , see [30] Recalling the definition of the symmetric graph Laplacian L s given by (2.8) and the relation between the spectral decomposition of W and the one of W in (2.17), we observe that a normalisation step now needs to be computed for obtaining the spectral decomposition of L s . Defining 1 L as the L-dimensional vector consisting of ones and 1 S−L analogously, we use (2.24) and start computing the nonnegative vector (2.25) Therefore, the matrices W XX and W XY can be normalised simply by considering: where the division is intended element-wise and ⊗ is the standard vector tensor product. A further step of normalisation is now needed since the approximated eigenvectors of W , i.e. the columns of the matrix V in (2.23) may not be orthogonal. Such normalisation may be obtained by using auxiliary unitary matrices. We refer the reader to [9, Section 3.2] for more details on this.
Once these additional normalisation steps are completed, we then get a spectral decomposition of W in terms of its eigenvaluesλ i and the corresponding normalised eigenvectors v i , i = 1, . . . , S. Therefore, recalling (2.17), the spectral decomposition of L s is given in terms of the eigenvalue 1 −λ i and eigenvectors v i . Algorithm 1 GL-minimisation with Nyström extension for image segmentation 1: Parameters: L S, σ, ε, C. 2: select L random points and build the set X ⊂ V 3: get a partition V = X ∪ Y, Y := X c 4: determine features and edge weights of X and Y using (2.10) and build W XX and W XY 5: Nyström extension to compute W XY (2.24) and matrix V of eigenvectors of W (2.23) 6: get normalisedŴ XX andŴ XY using (2.25)-(2.26) 7: normalise V and get eigenvalues-eigenvectors of W (λ i , v i ) 8: output ← eigenvalues-eigenvectors (1 −λ i , v i ) of L s used as GL minimisation input 9: convex splitting for GL minimisation through Fourier transform methods, as described in Section 2.3 10: output ← the binary segmentation.

Pseudocode
We present here the pseudocode combining all the different steps described above for the realisation of the GL minimisation. We recall that ε is the scale parameter of the GL functional (2.9), σ is the variance used in the Gaussian similarity function (2.10), C is the convex splitting parameter in (2.14a)-(2.14b) and L is the number of sample points in (2.19).
We will now give further details. First we randomly select L pixels from I. As described in Section 2.2 we now create a vertex set V ∼ = I, which we partition into a set X, consisting of the vertices corresponding to the L randomly chosen pixels, and a set Y := V \ X. We now compute the feature vectors of each vertex in V . If I is a grey scale image, we can represent features by an intensity map f : V → R. If I is an RGB colour image instead, we use a vector-valued (red, green, and blue) intensity map f : We mirror the boundary to define neighbourhoods also on the image edges. The feature function ψ : V → R K concatenates the intensity values in the neighbourhood of a pixel into a vector: ψ(x) := (f (ν 1 (x)), . . . , f (ντ (x))) T , where ν(x) = (ν 1 (x), . . . , ντ (x)) ∈ Rτ is the neighbourhood vector of x ∈ V andτ = (2τ + 1) 2 , the size of the neighbourhood of x. Note that K =τ if I is a grey scale image and K = 3τ if I is an RGB colour image.
Using (2.10), the matrices W XX and W XY from (2.22) can now be constructed and the Nyström extension can be performed for approximating the eigenvectors and eigenvalues of W as described in Section 2.4, which are then used to compute the eigenvectors {v k } of L s and corresponding eigenvalues {λ k }, compare (2.17). Recalling (2.16), those eigenvectors are used as basis functions for U n , the numerical approximation of u in the n-th iteration of the GL minimisation. Considering (2.15) and writing the nonlinear quantities appearing in terms of {v k } similarly as in (2.16), we have The computation of U in the next iteration reduces to finding the coefficients α k n+1 in the expression in terms of β k n , γ k n and the other parameters involved, that is the scale parameter ε in (2.9), the parameter C > 0 appearing in the splitting (2.14) and the time step ∆t. Using (2.15), we compute α k n+1 simply as where D k is defined as D k := 1 + ∆t(ελ k + C).

The Hough transform
In this section we give a brief description of a mathematical transformation that detects objects in an image with specific, a priori specified shapes. For our purposes, we focus on straight lines (for which the Hough transform was originally introduced and considered [37]) and circles, [26]. Other applications of this transformation for more general curves exist as well. In [7,42] the Hough transform is used in the context of astronomical and medical images for a specific class of curves (Lamet and Watt curves). In [34] applications to cellular mitosis are presented. There, the Hough transform recognises the cells (as circular/elliptical objects) and tracks them in the process of cellular division. The general idea behind the use of the Hough transform is to map the ambient space to an auxiliary space, called the parameter space (as it is related to the parametric representation of the geometrical objects we are interested in). There, objects with specified shapes are easily recognisable as peaks of specific functions. Let us clarify these concepts with two examples.
Detecting line segments. We start from the typical slope-intercept form of a line: y = mx + q, m, q, x, y ∈ R. (3.1) Traditionally, the equation above is considered as a function of points with coordinates (x, y) satisfying equation  Hence, if we are given a black and white image in the x-y plane, and for all coordinates (x, y) of black locations in the image, we draw the corresponding lines in the m-q plane, intersection points of those lines will tell us which (x, y) locations in the image lie on the same line. Of course any two points lie on a line, thus we are specifically interested in intersection points in the m-q plane in which many different lines intersect, indicating the presence of an actual black line segment in the original image.
Drawbacks of this parametrisation are the need for an unbounded parameter space to describe near vertical lines and the impossibility to describe a vertical line. One alternative is the normal parametrisation which views a straight line in x-y space as the tangent line to a circle with radius ρ, touching the circle at angular coordinate θ, as illustrated as in Figure 5a, [26]. In ρ-θ parameter space this leads to ρ = x cos θ + y sin θ, θ ∈ [0, π]. (3. 2) The objects in the parameter space are now sinusoidal curves, but again intersection points identify parameters for the points lying on the same straight line in the x-y plane. Figures 5b and 5c show a binary image with two black straight lines and the corresponding parameter space. The bright spots in the parameter space indicate a large number of intersections, thus identifying the two lines in the original image. With parametrisation (3.2) the Hough transform can be interpreted in terms of the Radon transform, as detailed in [54]. The Radon transform describes mathematically several medical imaging acquisition problems (PET and CT, for example) and, as such, is widely used in their mathematical modelling. We refer the reader to [53,16] and the references therein for a detailed variational analysis of such problems in appropriate measurement and image spaces.
Detecting circles. Analogously to what we did above, when looking for circular structures in a given image, we consider, for (x, y) ∈ R 2 the parametric representation of a circle,    Numerical strategy. We begin by using the Canny method (see A.1 in the Appendix) to identify the edges in the image. Having decided which geometrical shape we are interested in (and, as such, its general parametric representation), the corresponding parameter space is subdivided into accumulator arrays (cells) whose dimension depends on the dimension of the parameter space itself (2D in the case of straight lines, 3D in the case of circles). Each accumulator array groups a range of parameter values. The accumulator array is initialised to 0 and incremented every time an object in the parameter space passes through the cell. In this way, one looks for the peaks over the set of accumulator arrays as they indicate a high value of intersecting objects for a specific cell. In other words, they are indicators of potential objects having the specific geometrical shape we are interested in.

Pseudocode
Numerically, dealing with the Hough transform consists of looking for peaks of the accumulator arrays in the parameter space onto which the original image is mapped. We use the MATLAB routines hough, houghpeaks, and houghlines for straight lines detection and imfindcircles for circle detection. The accuracy and the number of detections for such routines can be tuned by some parameters, such as, for instance, the maximum number of peaks one wants to consider, obj max , or the array peak threshold, thresh, i.e. the minimum number of elements for an accumulator array to be considered a peak. The user also has to specify an initial range of pixel values [s min , s max ] as a very rough approximation of the measurement scale. Namely, in the case of line detection this determines a minimum/maximum spacing between lines , whereas for circle detection this serves as a rough approximation of the range of values for the circles' radii. This rough approximation may be given, for example, from average data which the user knows a priori. We explain this with some examples in Section 4. Accuracy of the detection algorithm is tuned by a parameter acc. In case of linear objects this corresponds to choose the maximum number of pixels between two line segments to consider them as one single line, whereas for circle detection this corresponds to the circularity of an object to be considered a circle. We report in this section the numerical results obtained by the combination of the methods presented for the detection and quantitative measurement of objects in an image.
To avoid confusion, we will distinguish in the following between two different meanings of scale. Namely, by image scale we denote the proportion between the real dimensions (length, width) of objects in the image and their corresponding dimensions quantified in pixel count. Dealing with measurement tools, we talk about measurement scale to intend the ratio between a fixed unit of measure (mm or cm) marked the measurement tool considered and the correspondent number of pixels on the image.

Male pied flycatcher's blaze segmentation
Here we present the numerical results obtained by applying algorithms 1 and 2 to the male pied flycatcher blaze segmentation and measurement problem described in the introduction. Our image database consists of 32 images of individuals from a particular population of flycatchers. Images are 3648 × 2736 pixels and have been taken by a Canon 350D camera with Canon zoom lens EFD 18-55 mm, see Figure 1. In each image one of two types of measurement tool is present: a standard ruler or a surface on which two concentric circles of fixed diameter (1 cm the inner one, 3 cm the outer one) are drawn. In the following we will refer to these tools as linear and circular ruler, respectively. Here, the measurement scale corresponds to the distance between ruler marks for linear rulers and to the radius of the inner circles for circular rulers. Figure 1 shows clearly that the scale of the images in the database may vary significantly because of the different positioning of the camera in front of the flycatcher. In order to study possible correlations between the dimensions (i.e. perimeter, area) of the blaze and significant behavioural factors, the task then is to segment the blaze and detect automatically the scale of the image considered to provide scale-independent measurements.
Comparison with Chan-Vese model. Due to the very irregular contours and the fine texture on the flycatcher's forehead, standard variational segmentation methods such as Canny edge detection (Section A.1) or Mumford-Shah models, [49,4,12,13], are not suitable for our task, as preliminary tests showed. Chan-Vese (Section A.2) is not suitable either, mainly because of the small scale detection limits, the dependence on the initial condition, and the parameter sensitivity which may prevent us from an automatic and accurate segmentation of the tiny, yet characteristic feathers composing the blaze. In particular, the optimal parameters µ and ν appearing in the Chan-Vese functional (A.2) and a sufficiently accurate initial condition need to be chosen typically by trial and error for every image at hand.
On the other hand, the graph GL method exploits similarities and differences between pixels in terms of the RGB intensities and texture of their neighbourhoods. In our image database these similarities and differences are very distinctive and will guide the segmentation step, with less dependence on the image scale and the initial condition. Recalling Section 2.5, we note that some parameters need to be tuned for the graph GL minimisation as well. They are the number of Nyström points, the variance σ of the similarity function (2.8), the GL parameter ε and the parameter C for the convex splitting (2.14). However, differently from what we discussed for the Chan-Vese model, in our numerical experiments we had to tune these parameters only once, observing a good level of robustness of the method for all the images in the database with respect to such a choice.
For comparison, we report in Figure 7 the blaze segmentation results obtained by using Chan-Vese model (see Section A.2 in the Appendix) and our graph-based method which will be described in more detail in the following.

Detailed description of the method
We divide our task into different steps: 1. For a given, unsegmented image, we detect the head of the pied flycatcher through a comparison with a user prepared dictionary (see Figure 8) using GL segmentation Algorithm 1. Further computations are restricted to the head only.
2. Starting from the reduced image, a second step similar to Step 1 is now performed for the segmentation of the blaze, using again Algorithm 1. A dictionary of blazes is used an extended set of features is considered.
3. A refinement step is now performed in order to reduce the outliers detected in the process of segmentation.
4. We use the Hough transform based Algorithm 2 to detect in the image objects with a priori known geometrical shape (lines for linear rulers, circles for circular rulers) for the computation of the measurement scale.
5. The final step is the measurement of the values we are interested in (i.e. the perimeter of the blaze, its area and the width and height of the different blaze components). In the case of linear rulers our results are given up to some error (due to the uncertainty in the detection of the measurement scale computed as average between ruler marks distances).
In order to establish relations with behavioural and biological data confirming or contradicting the initial assumption of correlation between blaze size and higher attractiveness presented in the introduction [52], we have implemented a user ready program for the quantitative analysis and measurements of the size of the bird blazes which is currently used by the Department of Zoology of the University of Cambridge. The results of this study will be the topic of a forthcoming paper [17].
In the following we give more details about each step.
Step 1: Head detection. We consider unlabelled images in the database and compare each of them with a dictionary of previously labelled images, see Figure 8. The training regions (i.e. the heads) are labelled with a value 1, the background with value −1. Unlabelled regions are initialised with value 0. The main computational difficulties in this step are due to size of the images considered. This may affect the performance of the algorithm as in order to apply the Nyström completion technique described in Section 2.4 one has to choose an adequate number of points whose features will approximate well the whole matrix. The larger and more heterogeneous the image is, the larger will be the number of points needed to produce a sensible approximation. We circumvent this issue noticing that at this stage of the algorithm, we only need a rough detection of the head which will be used in the following for the accurate segmentation step. Thus, downscaling the image to a lower resolution (in our practice, reducing the resolution by ten times the original one) allows us to use a small number of Nyström sample points (typically 150-200) to produce an accurate result.
For this first step we use as features simply the RGB intensities and proceed as described in Section 2.5. Once the head is detected, the resulting image is upscaled again to its original resolution. The solutions computed for the images in Figure 1 are presented in Figure 9.
Step 2: Blaze segmentation. We consider now the reduced image from which we want to extract the flycatcher's blaze. Again, a dictionary of different blazes is manually created by the user (see Figure 10). Again, training regions (the blazes) are labelled with value 1 and the black feathers in the background with value −1.
As before, unlabelled regions are initialised with value 0. At this stage, RGB intensities alone are not enough to differentiate the blazes from the background consistently in a large number of bird images, due to the colour difference between different blazes. For this step, an additional feature to be considered is the texture of the blaze. For this purpose, we use the MR8 texture features presented in [65] and proceed as detailed in Section 2.5. For 3×3 neighbourhoods, the feature vector for each pixel will be an element in R 99 , see Section 2.5. The Ginzburg-Landau minimisation provides the segmentation results shown in Figure 11.   Step 3: Segmentation refinement. This step uses very simple morphological operations in order to remove false detections obtained after Step 2. These can occur due to the choice of colour-texture based features used to compute the feature vectors in Step 2. For instance, when looking at Figure 11 (right) we observe that some bits on the left pied flycatcher's cheek have been detected as they exhibit similar texture properties as the ones on the blaze. In order to prevent this, our software asks the user to confirm whether the segmentation result provided is the expected one or if there are additional unwanted regions detected. If that is the case, using the MATLAB routine bwconncomp we label all the connected components segmented in the previous step, discarding among them all the ones whose area is smaller than a fixed percentage (typically 10%) of the largest detected component (presumably, the blaze). This works well in practice, see Figure 12. If the user is not satisfied he or she can remove manually the unwanted regions. Figure 13 shows some blaze segmentation results after the refinement step.  Step 4: Measurement scale detection. The images in our database divide into two groups: the first is characterised by the presence of linear rulers, whereas the second contains circular rulers (Figure 1). We thus need to use the Hough transform based Algorithm 2 to detect lines or circles, respectively. The user is then required to tell the software which objects he or she wants to detect. In both cases, in order to avoid false detections (such as "aligned" objects erroneously detected as lines, or circle-like objects wrongly considered as circles, see Figure  14), a good candidate for a rough, sensible approximation of the measurement scale is needed as described in Section 3.1. In order to get this, we proceed as follows: after detecting the head as in Step 1, we use the option EquivDiam of the MATLAB routine regionprops to detect the diameter of the head region (in pixels). We then compare such measurement with pre-collected average measurements of head diameters of male pied flycatchers of a similar population (in cm), thus obtaining an initial approximation of the measurement scale. In the case of images containing linear rulers, this will serve as a spacing parameter s for the algorithm. In other words, only lines distant at least s pixels from each other will be considered. In the case of circular rulers, the same rough approximation will serve similarly as an indication of the range of values in which the Hough transform based MATLAB function imfindcircles will look for circles' radii. For linear ruler images, the algorithm will look only for parallel lines aligned with a prescribed direction. We set this direction as the one perpendicular to the longest line in the image (since the expectation is, that this longest line is the edge of the ruler). Results of this step are shown in Figure 15. Outliers removal for linear rulers. The scale detection step described above may miss some lines on the ruler. This can be due to an oversmoothing in the denoising step, to high threshold values for edge detection or also to the choice of a large spacing parameter. Furthermore, as we can see from Figures 1 and 15, we can reasonably assume that the ruler lies on a plane, but its bending can distort some distances between lines. Moreover, few other false line detections can occur (like the number 11 marked on the ruler main body in Figure 15). To exclude these cases, we compute the distance (in pixels) between all the consecutive lines detected and eliminate possible outliers using the standard interquartile range (IQR) formula [62] for outliers' removal. Indicating by Q 1 and Q 3 the lower quartile and the third quartile, an outlier is every element not contained in the interval Finally, we compute the empirical mean, variance and standard Step 5: Measurement. Once the measurement scale has been detected, it is easy to get all the required measurements. We are interested in the perimeter, the area of the blaze and also in the height and width of the whole blaze component. For linear rulers, due to the error committed in the scale detection step, these values present some uncertainty and variability (see above). In Table 1 we show the results of numerical tests on a sample of 30 images with linear rulers. For every image in the sample we compute the standard deviation (SD) error and report in the table the minimum, maximum, and average SD error over the single ones compute, together with the relative standard deviation (RSD) which gives a percentage indication of the error committed.
where σ is the sample SD andX is the sample mean of measurements. We observe a minimum and maximum SD of 4.00 and 10.67 pixels, respectively, which, compared to the dimension of the original image (3648 × 2736 pixels) suggests a reasonable precision. This is confirmed by the average SD value over the sample which is found to be 6.81 pixels. In percentage, the average error over the sample is 11.99%. For circular rulers, we observed in all our experiments that an initial approximation of the range of values for the circle radius (see Step 4 above) results in a robust and typically outlier-free detection of the circular ruler and consequently in an accurate measurement of its radius; the only possible cause of variability and error is its bending. Uncertainty in the measurements of lengths and areas is calculated with standard formulas in propagation of errors.

SD min
SD max mean SD RSD min RSD max mean RSD 4.01 pixels 10.67 pixels 6.81 pixels 6.59 % 17.36 % 11.99 % Table 1: Precision of the measurement scale detection for linear rulers on a sample of 30 images. The minimum, maximum and average standard deviation (SD) error together with the corresponding relative standard deviation (RSD) errors are reported.
Despite these variabilities, our method is a flexible and semi-supervised approach for this type of problem. Further tests on the whole set of images and improvements on its accuracy are a matter of future research. The analysis of the resulting data measurements for the particular problem of flycatchers' blaze segmentation will be the topic of the forthcoming paper [17].
We compare in Table 2 between the use of our combined approach and the use of the manual line tool of the IMAGEJ software for the measurement of the blaze area. Namely, we measured in Figure 1b and in Figure 1c the ruler scale by means of the IMAGEJ line tool by considering, for each image, two different 3 cm-sections of the ruler; we then measured manually the number of pixels each section crossed, divided each measurement by 30 and averaged the two results to obtain an estimate of the ruler scale (i.e. the number of pixels crossed by a 1 mm horizontal or vertical line segment). We then measured the area of the blaze after segmenting it by means of the 'magic-wand' [48] IMAGEJ tool and trapezium fitting [52] (see Figure 2). The results are reported in Table 2 both as number of image pixels inside the blaze and in mm 2 , where this second value has been calculated using the measurement scale detected as described above. We then repeated such measurements using our fully automated Hough transform method for ruler scale detection, reporting as above the measurements of the blaze area computed both as number of image pixels and in mm 2 . We observe a good level of accuracy of our combined method (see also Table 1) with respect to the 'magic-wand' manual approach of Moreno [48], while, unsurprisingly, the blaze measurements obtained by pure trapezium fitting as proposed by Potti and Montalvo in [52] tend to overestimate the area of the blaze.  Table 2: Comparison between ruler scale detection by using manual IMAGEJ line tool and our Hough Transform (HT) method with corresponding measurements of the segmented blaze area obtained by using IMAGEJ 'magicwand' (MW) tool [48], trapezium fitting (Trap.) [52] (see also Figure 2) and the graph Ginzburg-Landau (GL) minimisation.

Moles monitoring for melanoma diagnosis and staging
In this section we focus on another application of the scale detection Algorithm 2 in the context of melanoma (skin cancer) monitoring, see Figure 3. Early signs of melanoma are sudden changes in existing moles and are encoded in the mnemonic ABCD rule. They are Asymmetry, irregular Borders, variegated Colour and Diameter ‡ . In the following we focus on the D sign. Due to their dimensions and their irregular shapes, moles are often very hard to measure. Typically, a common dermatological practice consists in positioning a ruler under the mole and then taking a picture with a professional camera. Sudden changes in the evolution of the mole are then observed by comparison between different pictures taken over time. Hence, their quantitative measurement may be an indication of a malignant evolution In the following examples reported in Figure 16, we use the graph segmentation approach described in algorithm 1 where texture-characteristic regions are present (see Figure 16a) and the Chan-Vese model (A.2) for images characterised by homogeneity of the mole and skin regions and the regularity of mole boundaries (Figures (16b)-(16c)). For the numerical implementation, we use the freely available online IPOL Chan-Vese segmentation code [31]. Let us point out here that previous works using variational models for accurate melanoma segmentation already exist in literature, see [20,1], but in those no measurement technique is considered.

Other applications: animal tracks and archeological finds' measurement
We conclude this section presenting some other applications for the combined segmentation and scale detection models presented above.
The first application is the identification and classification of animals living in a given area through their soil, snow and mud footprints. Their quantitative measurement is also interesting in the study of the age and size a of a particular animal species. As in the problems above, such measurement very often reduces to a very inaccurate measurement performed with a ruler placed next to the footprint image. In Figure 17a § our combined method is applied for the measurement of a white-tailed deer footprint.
As a final application, we focus on archaeology. In many archaeological finds, objects need to be measured for comparisons and historical studies [36]. Figure 17b shows the application of our method to coin measurements. Due to its circular shape, for this image a combined Hough transform method for circle and line detection has been used. The example image is taken from [36] where the authors propose a gradient threshold based method combined with a Fourier transform approach. Despite being quite efficient for the particular applications considered, such approach relies in practice on the good experimental setting in which the image is taken: almost noise-free images and very regular objects with sharp boundaries (mainly coins) and homogeneous backgrounds are considered. Furthermore, results are reported only for rulers with vertical orientation and no bending. ‡ Prevention: ABCD's of Melanoma. American Melanoma Foundation, http://www.melanomafoundation.org/prevention/abcd. htm.

Conclusions
In this paper we consider image segmentation applications involving measurement of a region's size, which has applications in several disciplines. For example, zoologists may be interested in quantitative measurements of some parts of the body of an animal, such as distinctive regions characterised by specific colours and texture, or in animal tracks to differentiate between individuals in the species. In medical applications, quantifying an evolving, possibly malignant, mass (like, for instance, skin melanoma) is crucial for an early diagnosis and treatment. In archaeology, finds need to be measured and classified. In all these applications, often a common measurement tool is juxtaposed to the region of interest and its measurement is simply read directly from the image. This practice is typically inaccurate and imprecise, due to the conditions in which pictures are taken. There may be noise corrupting the image, the object to be measured may be hard to distinguish, and the measurement tool can be misplaced and far from the object to measure. Moreover, the scale of the image depends on the image itself due to the varying distance from the camera of the ruler and objects to measure. The method presented (based on [9]) consists of a semi-supervised approach which, by training the algorithm with some examples provided by the user, extracts relevant features from the training image (such as RGB intensities, texture) and uses them to detect similar regions in the unknown image. Mathematically, this translates into the minimisation of the discrete Ginzburg-Landau functional defined on graphs. To overcome the computational issues due to the size of the data, Nyström matrix completion techniques are used and for the design of an efficient numerical scheme, convex splitting is applied. The measurement scale detection task is performed by using the Hough transform, a geometrical transformation which is capable of detecting objects with a priori known geometrical shapes (like lines on a ruler or circles with fixed diameter). Once the measurement scale is detected, all the measurements are converted into a unit of measure which is not image-dependent.
Our method represents a systematic and reliable combination of segmentation approaches applied to several real-world image quantification tasks. The use of dictionaries, moreover, allows for flexibility as, whenever needed, the training database can be updated. With respect to recent developments [68] in the fields of data mining for the analysis of big data, where predictions are often performed using training sets and clustering, our approach represents an interesting alternative to standard machine learning (such as k-means) algorithms.
3. Suppression: only local maxima are marked as edges. In this step the algorithm thin the blurred edges detected in Step 2 in order to get sharper edges. This translates into discarding all the points in the gradient image where there are not local maxima.
4. Double thresholding: this step refines the suppression Step 3 by using two thresholding values to identify strong edges (i.e. with gradient magnitude larger than the higher threshold), weak edges (gradient magnitude between the two thresholds) and to discard false edges (gradient magnitude lower than the lower threshold).

5.
Edge tracking by hysteresis: false edges not connected to strong edges from the previous step are discarded. Only the final, potential edges are maintained and the final edge set is determined.
Despite being very easy and widely applicable, the Canny algorithm described above has a serious drawback: many false detections can occur, for example if too much noise is left in the image after Step 1, or if texture is mistaken for edges.

A.2 The Chan-Vese model
To identify edges of an image, one can approximate the original image by a piecewise regular function whose changes in regularity correspond with the edges of the image. In 1989 Mumford and Shah were the first to propose a variational method that approximates the original image f with a piecewise smooth function u. Due to the lack of differentiability of the model, smoother approaches have been considered in [4,12,13] where also Γ-limits are studied.
Another well-known model is the Chan-Vese model presented in [21] which can be considered as a simplification of the Mumford-Shah model where the image u we seek is assumed a priori to be binary: In the definition above, C is the boundary of a closed set. The variational Chan-Vese model seeks the optimal c 1 and c 2 (i.e. the optimal u of the form (A.1)) and contour C to minimise F(c 1 , c 2 , C) := µ Length(C) + ν Area(int(C)) + λ 1 where µ, ν, λ 1 , λ 2 > 0. In (A.2) the first term penalises the length of C, ensuring its regularity. The size of C is controlled by the second term which is penalising the area in the interior of C, while the two other terms penalise the discrepancy between the piecewise constant u and the given image f in the interior and exterior of C, respectively. By computing the local minima of (A.2) one retrieves the optimal binary approximations u of f . Existence and regularity of solutions of (A.2) can be deduced from similar results for the solution of a simplified Mumford-Shah model [47]. If C is the zero level set of a Lipschitz continuous function φ(x), then the energy function F in (A.2) can be rewritten as function of φ as µ Ω δ(φ(x))|∇φ(x)| dx + ν Ω H(φ(x)) dx+ where H denotes the Heaviside function and δ its distributional derivative, i.e. the Dirac mass. For numerical implementations, H and δ are usually regularised. In order to find the minima of (A.3) an alternating strategy is used, updating c 1 , c 2 and φ successively in each iteration. The level set representation of C is initially constructed as Thus, the inside and the outside of C are distinguished by the sign of φ. Several extensions of the Chan-Vese model have been considered in literature. In particular, in [22] the model is extended for vector-valued data, the only difference with (A.2) being the use of Euclidean vector norms in the integrands.
Despite being very popular and widely used in applications, the Chan-Vese model and its extensions present intrinsic limitations. Firstly, due to the lack of uniqueness of solutions of the minimisation problem, the segmentation result is strongly dependent on the initial condition. This is a significant drawback for automation of the