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