Modeling sets of unordered points using highly eccentric ellipses

Abstract

An algorithm for modeling a set of unordered two-dimensional points by line segments is presented. The points are modeled by highly eccentric ellipses, and line segments are extracted by the major axes of these elongated ellipses. At first, a single ellipse is fitted to points which is then iteratively split to a large number of highly eccentric ellipses to cover the set of points. Then, a merge process follows in order to combine neighboring ellipses with almost collinear major axes to reduce the complexity of the model. Experimental results on various databases show that the proposed scheme is an efficient technique for modeling unordered sets of points and shapes by line segments. A computer vision application of the method is also presented regarding the detection of retinal fundus image features, such as end-points, junctions, and crossovers.

1 Introduction

In many computer vision applications, at a mid-level process, it is common to fit line segments in order to model a set of unordered points so as to summarize higher level features. For example, the detection of vanishing points [1], the vectorization of raster images [2], and the detection of road structures and parts [3] are among applications necessitating line segment description of image structures. In many of the aforementioned problems, the involved algorithms assume that they are provided with an ordered point set and standard polygonal approximation [4, 5] is then applied. However, determining the ordering of point sets is not a trivial task, and in the method described herein, we relax this assumption by making no prior hypothesis about the ordering of the points.

In the above context, the Hough transform (HT) is a widely used method for line fitting, and many variants have been proposed to improve its efficiency [6, 7]. One of these variants is the randomized Hough transform (RHT) [8, 9] which randomly selects a number of pixels from the input image and maps them into one point in the parameter space which was shown to be less complex, compared to the original algorithm, as far as time and storage issues are concerned. In [10], the probabilistic HT was proposed whose basic idea is to apply a random sampling of edge points to reduce computational complexity and execution time. Further improvements were introduced in [11]. A similar concept was proposed in [12], where an orientation-based strategy was adopted to filter out inappropriate edge pixels, before performing the standard HT line detection which improves the randomized detection process. Also, the idea of fuzziness is integrated in the main algorithm in [13] to model the uncertainty imposed to the contour due to noise. Thus, a point can contribute to more than one bin in the standard HT process. A general comparison between probabilistic and non-probabilistic HT variants can be found in [14].

The robust HT is introduced in [15] where both the length and the end-points of the lines may be computed. Moreover, the algorithm in [16] provides a method for adopting a shape-dependent voting scheme for the calculation of the histogram bins. Finally, a novel HT based on the eliminating particle swarm optimization (EPSO) is proposed in [17], to improve the execution time of the algorithm. The problem parameters are considered to be the particle positions and the EPSO algorithm searches the optimum solution by eliminating the ‘weakest’ particles, to speed up the process.

Line segment fitting may also be used in a shape description process. The commonly used algorithm of Moore [18] was a first solution to shape following. However, this algorithm is appropriate only for traversing curves without intersections and produces models with high complexity, although improvements of the main algorithm have also been considered up to date [19]. Another common model fitting method is the RANSAC algorithm [20], which despite the fact that it provides robust estimations, it is appropriate for fitting only one model at a time. Other approaches are the incremental line fitting [21] which is sensitive to noise and, most importantly, needs sequential ordering of the points and probabilistic methods [22] based on the Expectation-Maximization algorithm, generally necessitating the prior determination of the number of model components.

In this paper, we propose a method to model a set of unordered points by line segments. The result of the method is a sequence of straight line segments modeling the points, which correspond to the major axes of highly eccentric ellipses fitted on them. We call the procedure Direct Split-and-Merge algorithm or DSaM. The points may correspond to the locations of edge pixels. The basic idea relies on a split-and-merge process, where the algorithm is initialized by one ellipse, representing the mean and the covariance of the initial point set, which is then split, through an iterative scheme, until a number of small and elongated ellipses occurs. Then, a merge process takes place to combine the resulting ellipses in order to reduce the complexity of the representation. At the end, the algorithm provides a relatively small number of elongated ellipses fitted on the points. The proposed method is general and may also model a set of points containing inner structures, such as joints. Numerical experiments on common shape databases are presented to underpin the performance of the proposed method. Also, edgemaps of natural images were used in our experiments to verify the efficiency of the method on real data. The number of ellipses depends on the structure of points, and it is determined by two parameters controlling both the split and the merge processes. The extension of the method to a three-dimensional (3D) case is trivial, as the elongated ellipses may be transformed to 3D ellipsoids. In the 3D case, our rationale is to model the points by highly eccentric ellipsoids, thus, instead of line segments, multiple 2D planes fit to points. In that case, data from range images could be modeled. A preliminary version of this method was introduced in [23]. Note that our method is not affected by degradations in the contour, and thus, it is more reliable compared to the distance transform [24] or the medial axis transformation [25], as it assumes more relaxed initial conditions (a set of unordered points) compared to conventional shape description methods, where point traversal should be known a priori.

Moreover, we present some applications that could benefit by our DSaM method. To that end, an efficient algorithm for retinal image analysis based on the proposed modeling framework is presented. More precisely, the purpose is to parse the image representing the retinal vessels and extract their characteristic points (end-points, junctions, and crossovers).

The remainder of the paper is organized as follows: In Section 2, the split-and-merge algorithm is presented. Experimental results using various databases (MPEG7 [26], GatorBait100 [27], ETHZ [28], Brown [29]) that contain either object silhouettes or collections of real images are shown in Section 3. In the same section, a comparison of the proposed method with a variant of the Hough transform, namely, the progressive probabilistic Hough transform (PPHT) is demonstrated [30]. In Section 4, we present an applications of our method: an algorithm for retinal image analysis is presented along with its evaluation using the DRIVE database [31–33].

2 Direct split-and-merge method

Let X = {x _i |i = 1,…,N} be a set of points and E = {ε _j |j = 1,…,K} be the set of line segments modeling the points, where ε _i is the line describing the i th segment.

We define the modeling error Δ induced by the representation of line segments:

$$\Delta (X,E)=\sum _{i=1}^N\sum _{j=1}^K{\delta }_{\mathit{\text{ij}}}d({\mathbf{x}}_i,{\epsilon }_j),$$

(1)

where K is the number of line segments the model uses to model the points, ${\mathbf{x}}_i\in {\mathbb{ℝ}}^2$, i = 1,…,N are the points, d(x _i ,ε _j ) is the perpendicular distance of point x _i to line ε _j , δ _{i j} is an indicator function whose value is one if point x _i belongs to line segment ε _j and is zero otherwise.

In order to prevent overfitting, models having a large number of line segments should be penalized. Therefore, an optimal model would have both low value of Δ and low complexity.

The computation of the ellipses, modeling the line segments, is performed in two steps: an iterative split process, where points are modeled by a number of line segments represented by the major axes of the corresponding ellipses and an iterative merge process, where small line segments are merged to reduce the model complexity. The split process tries to minimize the modeling error while the merge process decreases the model complexity, i.e., the number of line segments compared to the total number of points in the set.

What follows are the two steps presented in detail.

2.1 Split process

The ultimate goal of this step is to cover the point space with line segments representing the long axes of elongated ellipses and therefore, each point of the shape should be assigned to an eccentric ellipse. A split criterion is defined, based on Gestalt theory [34], which models the linearity and the connectivity the human brain uses when modeling contours.

In order to split a set X, it should be either non-linear or disconnected, or both. Linearity describes how close the points are to a straight line, while disconnectivity measures how concentrated these points are. In the ideal case, the covariance matrix of collinear 2D points should have a very large eigenvalue and a zero eigenvalue. The eigenvector corresponding to the larger eigenvalue indicates the direction of the line segment. If the linearity property is relaxed, the less collinear the points become (i.e., they diverge from the linear assumption) the larger the value of the minimum eigenvalue is. Based on that observation, in our method, linearity is described by the minimum eigenvalue of the covariance matrix of the points in X. Also, the disconnectivity W of two sets of points X, Y is the smallest distance between a point in X and a point in Y:

$$W(X,Y)=\underset{\begin{array}{c}\mathbf{x}\in X\\ \mathbf{y}\in Y\end{array}}{\min }|\mathbf{x}-\mathbf{y}|.$$

(2)

In the case of a single set, disconnectivity is the largest distance between two successive points in that set. It may be computed by projecting the points onto both axes defined by the eigenvectors of the covariance matrix of the set. Then, successive points are defined by scanning along the axes and their distances are computed. Let X _i be the projection of a set X onto the the eigenvector e _i . The disconnectivity of X is defined as

$$W\left(X\right)=\underset{\begin{array}{c}j=1,\dots ,N-1\\ i=1,\dots ,d\end{array}}{\max }|{\mathbf{x}}_i^j-{\mathbf{x}}_i^{j+1}|,$$

(3)

where N is the number of points in X, d is the dimension of X (here d = 2), and ${\mathbf{x}}_i^j$ is the j th point of the sorted set X _i . A large value of disconnectivity indicates a better separation of the point sets. The projections onto all of the eigenvectors should be examined as we do not know a priori which direction to follow while splitting. Although intuitively one would suggest to split along the direction of the principal axis, we observed that in many cases that approach was not the best. Also, let us note that as the ordering of the points is not known a priori, their projection onto the eigenvectors of their covariance matrix, provides a natural way of ordering.

The disconectivity of a single set of points is also important to be estimated in the split step, as there may exist subsets that although they are linear, they are disconnected, as it is demonstrated in Figure 1.

Figure 1

An example of an instance of a split where a subset of the points is linear but disconnected. If the disconnectivity of the points was not considered a single line segment would summarize these points.

The split of an ellipse should be performed along the direction defined by an eigenvector of its covariance. In order to select the split direction, the axis corresponding to an eigenvector is considered as the discrimination border between the split line clusters and points belonging to the same subplane are grouped together. Then, the disconnectivity of each line cluster is computed. Finally, the direction with the largest disconnectivity is selected for splitting (Figure 2).

Figure 2

Split process. (a) At iteration t + 1, the ellipse with center μ^t is split into two new ellipses e₁ and e₂, with centers ${\mu }_1^{t+1}$ and ${\mu }_2^{t+1}$ given by (4). (b) The new centers are marked with a star (*). The reassignment of the points to the new centers is shown. Points of one category, assigned to e₁, are marked with a square, while points assigned to e₂ are marked with a circle.

Eventually, the adopted strategy that minimizes Δ and prefers elongated ellipses can be expressed as follows: split every ellipse whose minimum eigenvalue is greater than a threshold T₁ (linearity) and the maximum gap, within the current segment is greater than a threshold T₂ (disconnectivity). The process is initialized with one ellipse, corresponding to the covariance of the initial points set centered at the mean value of the point locations. Thresholds T₁ and T₂ may be computed with a heuristic algorithm, as explained in Section 3.1.

At iteration t + 1, a given ellipse, characterized by the eigenvalues ${\lambda }_1^t$ and ${\lambda }_2^t$ of its covariance matrix Σ^t (with ${\lambda }_1^t\ge {\lambda }_2^t$), with center μ^t, is split to two new ellipses with centers the two antipodal points on the major axis:

$$\begin{array}{cc}{\mu }_1^{t+1}& ={\mu }^t+\sqrt{{\lambda }^t}{e}^t,\\ {\mu }_2^{t+1}& ={\mu }^t-\sqrt{{\lambda }^t}{e}^t,\end{array}$$

(4)

where e^t, λ^t are the eigenvector and the eigenvalue corresponding to the split direction along which split is performed (Figure 2).

The points of the split ellipse are then reassigned to the two new ellipses according to the nearest neighbor rule. In this way, new ellipses occur, which are more elongated as they have greater eccentricity and their minor axes are closer to the contour (Figure 3). Moreover, this detailed representation of the point set provides accurate modeling of the joints, corners, and parts of the contour exhibiting high curvature.

Figure 3

Steps of the split-and-merge process. The process is initialized with the mean and the covariance of the full set of points. (a) Split into two ellipses. (b) Split into four ellipses. (c) End of split (35 ellipses). (d) The final merge result (23 ellipses). The figure is better seen in color.

A variant of the method would be to compute the covariance matrix of the points on the convex hull of the point set, which provide more robustness to outliers.

2.2 Merge process

The role of the merge process is to reduce the complexity of the model. In case there exist adjacent ellipses whose major axes have similar orientations, it would be beneficial to merge and replace them with a more elongated ellipse. Therefore, in this step, ellipses are merged using the following rule: merge two consecutive ellipses, if the resulting ellipse has minimum eigenvalue smaller than a threshold T₁ (linearity) and the marginal width between the two line clusters is smaller than a threshold T₂ (disconnectivity).

Note that the threshold T₁ could be set equal to the threshold used in the split process, where the value of parameter T₁ specifies whether an ellipse has low eccentricity and needs to be split. In the merge process, it indicates whether two candidates for merging ellipses would result in an ellipse with high eccentricity. One could use the same threshold in both processes, assuming the same significance. On the other hand, a relaxation of the merge threshold could lead to a rougher model of the points, smoothing out details like joints. In our experiments, the merge threshold was selected to be the same with the split threshold. The same applies for threshold T₂ that indicates whether two segments are close enough to be considered as one line segment.

The overall description of the method is presented in Algorithm 1.

Algorithm 1 Direct Split-and-Merge Algorithm

3 Experimental results and discussion

In this section we evaluate the efficiency of the introduced algorithm. To that end, two categories of experiments were conducted. The purpose was to investigate the performance of the method both in shape data, but also in real images. Thus, various well-known databases were employed, that contain either object silhouettes or scenes of real images. The GatorBait100 database [27] consists of 38 shapes of different fishes grouped in 8 categories. The shapes of this database are not closed and contain many junctions. The MPEG7 shape database [26] consists of 1400 silhouettes of various objects clustered in 70 categories.The shape silhouette database used in [29], that contains 137 silhouettes of various objects, clustered in 13 categories, was also used in our experiments. Finally, to investigate the behavior of the proposed algorithm in real scene images, the images (257) from the ETHZ image set [28] were also used. Table 1 gives a brief description of each database. In all cases, the edges were extracted and the coordinates of the edge pixels were used to describe the contour. The Canny edge detector [35] was used in all cases.

Table 1 Short description of the databases used in our experiments

3.1 Numerical evaluation

In this section, we present the results^a of comparing the DSaM method with the widely used implementation of Kovesi [36]. Tables 2 and 3 summarize the numerical results, while Table 4 demonstrates the execution time for the experiments on MPEG7 [26], GatorBait100 [27], Brown [29], and ETHZ [28] datasets. Some representative images from those databases are given in Figure 4. As it can be observed, in some cases, there exist inner structures and thus, the ordering of the points is not obvious. Note that to share common parameters, in the Kovesi [36] implementation, we used the disconnectivity threshold of our method. The execution time for computing that value is not included in the execution time of the Kovesi implementation. The model complexity is computed by the index:

Table 2 Modeling error Δ ( 1 )

Table 3 Model complexity MC ( 5 )

Table 4 Average execution time (in seconds)

Figure 4

Some representative images of the databases we used in our experiments. Please note that in some cases inner structures exist. This does not permit to extract an ordering of the points (a) MPEG7 [26], (b) Gatorbait [27], (c) Brown [29], and (d) ETHZ [28].

$$\text{MC}=\frac{\#\text{ellipses}}{\#\text{points}}.$$

(5)

Lower values of MC imply lower complexity and therefore a more compact representation.

The distortion is the measure of the quality of the fitting and is computed as the average distance between a point and its corresponding line segment, as computed by each method. Please note that the average length of the diagonal of the bounding box of the various datasets is about 500 units (ranging from 300 units in Brown [29] to 700 units in ETHZ [28]).

In the proposed algorithm, there are two parameters to be a priori specified, a threshold that determines the elongation of an ellipse (T₁) and a threshold characterizing the disconnectivity of a set (T₂). Both parameters are used to decide whether to split (in the split process) or merge (in the merge process). A small value preserves the details, while a larger one provides more coarse results. For our experiments, we computed the values of the parameters as:

$$T_1=\frac{1}{N}\sum _{i=1}^N{\lambda }_i$$

(6)

$$T_2=\frac{1}{N}\sum _{i=1}^N{d}_i,$$

(7)

where N is the number of the points of the set, λ _i is the smallest eigenvalue of the covariance matrix of points $\left\{\mathbf{x}\right|\mathbf{x}\in N_{{\mathbf{x}}_i}^{\alpha },i=1,\dots ,N\}$, with $N_{\mathbf{x}}^{\alpha }$ being the α- neighborhood of x, and

$$d_i=\underset{\mathbf{y}\in N_{{\mathbf{x}}_i}^{\alpha }}{\min }\left|\right|{\mathbf{x}}_i-\mathbf{y}\left|\right|,i=1,\dots ,\mathrm{N.}$$

(8)

Large values for α decrease the model complexity providing larger modeling error and details are not preserved. In our experiments, we set α = ⌈0.01 × N⌉ for computing the values of T₁ and T₂. To make our implementation more efficient, instead of taking all points into consideration, we computed a random permutation of the indices of points and used only the first 10% of them. Thus, in high density datasets, like in the ETHZ database [28], the values of the thresholds could be estimated quickly.

In general, the DSaM method and the Kovesi implementation produce models with similar complexity, a fact that is obvious, since they employ the same thresholds. However, the DSaM method provides much more accurate results w.r.t distortion (Tables 2 and 3). Concerning the execution time of the methods, one may observe that for simple datasets, like the object silhouettes [26, 29], both methods are quite fast. However, for more complex datasets, like the ETHZ [28], DSaM is much faster (up to about eight times) than the Kovesi implementation, since these datasets, contain a lot of points and Kovesi implementation has to check them all, one by one.

As our method models the line segments with ellipses, we tried to fit line segments by exploiting various modifications of a typical Gaussian Mixture Model (GMM) [37], for example, using an incremental GMM, or imposing constraints in the update step of the covariance matrices (decomposing the covariance matrices with SVD, replacing the corresponding minimum eigenvalue with a very small value, threshold T₁, and then recomputing the covariance matrix). All these variants failed to produce an efficient result. Moreover, the execution time was quite high (ten times compared to those of DSaM). Thus, we opted for excluding the results from this presentation.

Finally, we conducted experiments to verify the robustness of the method against the presence of noise that degrades the contour of an object. To that end, we used three patterns (see Figure 5a,b,c) which were randomly repeated to create new images. A representative image is given in Figure 6. As a ground truth for the number of line segments, we used 4 for the square, 3 for triangle, and 10 for the star.

Figure 5

Images used to create the artificial dataset and a representative result of a degraded contour. (a, b, c) The primitive images used to create the artificial dataset for experiments with Gaussian additive noise. (d) Contour degraded by additive Gaussian noise of 18 dB.

Figure 6

A representative test image produced by randomly repeating the patterns of Figure 5 a,b,c.

The set of unordered points was produced by simple edge detection. Then, zero mean Gaussian noise with varying standard deviation was added in order to get several configurations of signal-to-noise ratio (SNR). A representative result of a degraded contour is given in Figure 5d. Note that no ordering of points may be established in that case and thus polygon approximation may not be performed. To make the experiment independent from the noise configuration, each experiment was repeated 20 times. The algorithm assumes that a form of binary data (e.g., an edge map) is provided. Degradation by noise is performed after the edge extraction in order to examine the behavior of the algorithm to the detection of line segments. If the noise was added to the original image, the edges would be erroneous and we would not have a standard baseline for evaluating the algorithm.

In Figure 7, we present the results of the experiments. The error is expressed as the absolute difference between the real number of segments and the one computed by our method. It can be observed that while the magnitude of the noise decreases, the error is also decreased. The difference between true and estimated number of segments is generally small, 3 on average with low variance (±2 segments), compared to the total number of line segments, 90 line segments on average, corresponding to 3% deviation between true and estimated measurement. Thus, it could be claimed that the proposed method exhibits a consistent and efficient performance even if the contour is corrupted.

Figure 7

Experimental results using the datasets of Figure6. They demonstrate the performance of our method in presence of Gaussian additive noise in terms of model complexity error. The vertical axis represents the absolute error between the real number of segments and the one computed by our method.

3.2 Comparison with the Hough transform

Since the proposed algorithm fits line segments to a set of points, we also tested it against the commonly used Hough Transform (HT). However, since the standard HT is appropriate for fitting lines and not line segments, we applied the Progressive Probabilistic Hough Transform (PPHT), as proposed in [30] and implemented in the OpenCV library [38]. The implementation of PPHT imposes three parameters: (1) a threshold, indicating the minimum number of points in a bin, in the line parameter space, in order to consider that the line is represented by a sufficient number of points, (2) the minimum length of a line segment, and (3) the maximum gap between line segments lying on the same axis. In our experiments, we fixed the last two parameters (after a trial and error procedure keeping those parameters that best fit the examined points) and varied the threshold. The obtained results for the PPHT exhibited significant irregularities such as a large number of overlapping lines for the same segment. Also, the corners of the shapes were not correctly captured. Representative experiments on the MPEG7 dataset [26] are shown in Figure 8a,b,c while the solution of our DSaM algorithm is illustrated in Figure 8(d).

Figure 8

Representative experiments on the MPEG7 dataset and the solution of our DSaM. (a,b,c) Results of the PPHT algorithm to a set of points representing the shape of a bone (MPEG7 dataset) y varying the minimum number of points in a bin (namely, 5, 15, and 25). Only a small fraction of the lines is drawn for visualization purposes. Note the overlapping lines. (d) The result of our method. The figure is better seen in color.

The PPHT is based on a histogram which correlates the accuracy of the result with the number of bins used. Also, a threshold must be established to eliminate lines with small participation in the final result. A small number of bins may lead to an underestimation of the number of segments, while a large number of bins increases the complexity of the model. As far as the threshold is concerned, its value may have similar effects in the final model. A large value may drop some segments, while a small value may be responsible for a large number of lines fitted, analogous to a GMM with one component per point. A more important drawback of the PPHT is that many overlapping lines may model the same line segment. Figure 8 presents solutions of PPHT for a given set of points and various parameters values.

4 Retinal fundus image feature characterization

In this section, we introduce an application of the proposed method. Note however that this is only an indicative presentation, as one could also use our algorithm for solving other problems that need the modeling of line segments, such as the extraction of road layers from raster maps [39], where the standard Hough Transform is used to detect lines in some part of their algorithm. Another application is the detection of the vanishing point in an indoor scene [40].

A methodology that extracts features from the retinal fundus image and characterizes them will be presented. The goal is to detect the intersection points between the vessels, as they could provide useful information to an automatic diagnosis system.

The detection and characterization of various topological features of the retinal vessels are an important step in retinal image processing algorithms within an autonomous diagnosis system. A deviation from common topological feature patterns may be an indicator of anomaly. A comprehensive study may be found in [41]. In a typical retinal vessel structure, more than 100 junctions may be present [32], a fact that makes the manual characterization a tedious and time-consuming task. Typical retinal features are presented in Figure 9.

Figure 9

The different features that the proposed algorithm can detect. The yellow point is an end-point, the orange point is an interior-point, and the green point is a crossover. All the other points are junctions (a T-junction is shown in red, and a bifurcation is shown in blue). The image is better viewed in color.

For our algorithm, three types of features are detected: end-points (points at the extremities of the vessels), interior-points (points along a vessel), junctions (a new vessel is a branch of a longer one - T-junctions or a vessel is split into two or more new vessels - bifurcation) and crossovers (one vessel passes over another). Please note that further processing is needed to distinguish between a crossover and a bifurcation.

To investigate the accuracy of the proposed algorithm, experiments were conducted on the DRIVE database [31], which includes 40 retinal images along with their manual extraction of the vessels. The ground truth used in [32], [33] was also employed. In our experiments, the manual segmentations were employed, as the scope of our algorithm is to detect junctions, crossovers, and end-points. The reader should refer to [42] or [43] for a detailed vessel extraction algorithm, which is a preprocessing step of the whole chain. At first, a Canny edge detector [35] is applied to extract the borders of the vessels and then a thinning algorithm [44] is used, to extract the center line of the vessels. In Figure 10a, the original image is shown. Figure 10b presents the manual segmentation, while the data used in our retinal parsing algorithm are shown in Figure 10c.

Figure 10

The original image, manual segmentation, data used in our retinal parsing algorithm, and confidence regions. (a) The original retinal image. (b) The manual segmentation of the vessels in (a). (c) Result of thinning the image in (b). (d) The confidence regions depicted as circles with a radius equal to 1% of the diagonal of the bounding box of the original set. The figure is better seen in color.

Once the center line of each vessel is extracted, the method presented in Section 2 is applied in order to extract the corresponding line segments and assign each point to its corresponding line segment. For each line segment, its extreme points are the points that have the largest distance from the corresponding extreme points of the same line cluster.

In Figure 11, the points x (summarizing the vessel structure) are depicted with red and black color, while the corresponding extreme points y are presented with green and blue dots. A rule is defined to characterize a point as end-point or junction or crossover. Let $\mathcal{C}\left(\mathbf{x}\right)$ be the index of the line cluster point where x belongs to. In Figure 11, two line clusters are shown ($\mathcal{C}\left(\mathbf{x}\right)=1$ and $\mathcal{C}\left(\mathbf{x}\right)=2$). To define the neighborhood of extreme points, a neighborhood radius threshold is defined as

$$T_n=w\ast \bar{d},$$

(9)

Figure 11

An instance of the point characterization algorithm. Points x (in red and black) correspond to the thinned lines of the extracted vessels. Green and blue points are the extreme points y computed by our DSaM algorithm. The yellow circle demonstrates the neighborhood of that point ($\mathcal{N}\left(\mathbf{y}\right)$). Red and black points lying in that circle are considered as neighbors of that extreme point. In that case, those points belong either to line cluster with index 1 or to line cluster with index 2. Thus,$\mathcal{C}\mathcal{N}\left(\mathbf{y}\right)=2$. The orange point corresponds to the nearest neighbor of y among the points of neighbor line cluster (black points). d is the minimum distance between the aforementioned nearest neighbor and the extreme points of its line cluster.

where $\bar{d}$ is the mean distance between all the nearest neighbors and w is a constant that is learned, as explained later in this section. Thus, for an extreme point y, the neighborhood $\mathcal{N}\left(\mathbf{y}\right)$ is the set of the points x such that ||x – y|| ≤ T _n . Note that $\mathbf{y}\in \mathcal{N}\left(\mathbf{y}\right)$. In order to characterize point y, we define $\mathcal{C}\mathcal{N}\left(\mathbf{y}\right)$ as the number of distinct line clusters that points $\mathbf{x}\in \mathcal{N}\left(\mathbf{y}\right)$ belong to.

In the example shown in Figure 11, the studied extreme point y is the green one, while its neighborhood $\mathcal{N}\left(\mathbf{y}\right)$ is defined by a circle centered at y with radius T _n . Red and black points lying within that circle are the neighbors of y. Those points belong either to line cluster 1 or to line cluster 2 and thus $\mathcal{C}\mathcal{N}\left(\mathbf{y}\right)=2$. If all neighbors of y belong to the same line cluster with that of y, then y would be an end-point ($\mathcal{C}\mathcal{N}\left(\mathbf{y}\right)=1$). In case where $\mathcal{C}\mathcal{N}\left(\mathbf{y}\right)>2$, y would be a junction or a crossover. The algorithm in its current form does not discriminate between them.

A special case occurs when $\mathcal{C}\mathcal{N}\left(\mathbf{y}\right)=2$, where the studied point y is either an interior-point or a junction (T-junction). In that case, further elaboration is needed to characterize the extreme point by examining whether $\mathcal{N}\left(\mathbf{y}\right)$ contains an extreme point or not. Thus, if y belongs also to the neighborhood of y^′, with y^′ denoting the neighbor of y, then y is an interior-point, otherwise it is a junction (T-junction). Please note that the neighboring relationship we are describing in that section is not reflective. For example in Figure 11, the green point, which is one of the extreme points of cluster 1, is neighbor to cluster 2 (as there are some points of cluster 2 within the yellow circle that defines the neighborhood of the green point). However, none of the extreme points of cluster 2 contain a point from cluster 1 in their neighborhood, and thus, cluster 1 is not neighboring to cluster 2.

In our example in Figure 11, this means that one of the two blue points would be inside the yellow circle. In case that $\mathcal{N}\left(\mathbf{y}\right)$ contains no extreme point, as shown in Figure 11, then we compute the minimum distance (denoted by d in Figure 11) between the nearest neighbor (orange point in Figure 11) of y among the points of the other neighboring line cluster (black points in Figure 11) and the corresponding extreme points (blue points in Figure 11). If $d\le \bar{d}$, then y is a (junction) (T-junction), otherwise it is an interior-point.

A detailed description of the rules used to characterize an extreme point y is presented in Algorithm 2.

Algorithm 2 Rules for vessel features characterization

Note that since the ground truth refers to the original vessels and not to their center lines, which is the input of our method, we determined a value T_conf that defines a confidence region around a ground truth point. A computed point is considered to match a ground truth point if it lies in its confidence region. In our experiments, T_conf is defined as a percentage (1%) of the length of the diagonal of the bounding box of points x. Figure 10d shows the confidence regions depicted as circles with a radius equal to T_conf. To establish a robust value for constant w ((9)), the precision and recall rates were computed for values of w between 1.8 and 4.0 with a step of 0.8. Then, the F measure (harmonic mean) was calculated as

$$F=2\frac{\mathit{\text{PR}}}{R+P},$$

(10)

where P is the precision and R is the recall:

$$P=\frac{\mathit{\text{TP}}}{\mathit{\text{TP}}+\mathit{\text{FP}}},$$

(11)

$$R=\frac{\mathit{\text{TP}}}{\mathit{\text{TP}}+\mathit{\text{FN}}},$$

(12)

where TP is the number of true positives, that is, the number of relevant items retrieved, FP is the number of false positives, that is, the number of irrelevant items retrieved and FN is the number of false negatives, that is, the number of relevant items not retrieved.

Figure 12 shows the plot of F measure for various values of parameter w in Equation (9). In our experiments, the F measure takes its maximum value for w = 2.9. The corresponding point is depicted with a black square in Figure 12. In that case, the corresponding precision is equal to 91.59% while the recall is 98.58%. The value of $\bar{d}$ computed from our experimental data is approximately equal to $\sqrt{2}$, which leads to T _n  = 4.10, Equation (9), corresponding to a neighborhood radius equal to 4 pixels.

Figure 12

The F measure, Equation ( 10 ), for various values of parameter w in Equation ( 9 ). The black square indicates the point that corresponds to the maximum value of F measure. This value (F=0.95) occurs for w=2.9 and provides a precision rate of 91.59% and a recall rate of 98.58%. More details are given in Section 4.

The mean execution time was 109 s for the extraction of the line segments and 12 s for the extraction and characterization of features using Matlab on a typical Dual Core 2 × 2.50 GHz machine with 2.0 GB RAM.

5 Conclusion

The problem of fitting several lines in order to describe a set of 2D unordered points is crucial in several computer vision applications. To that end, a split-and-merge algorithm was presented in this paper that iteratively fits a number of elongated ellipses to 2D points (DSaM). The points are represented by the major axes of the ellipses, while the minor axes account for possible perturbations from the linear modeling. Furthermore, a merge process combines neighboring ellipses with collinear major axes in order to reduce their number and consequently the complexity of the solution. An important aspect of the algorithm is that it does not require an ordering of the points and accurately handles contours or edge chains that contain joints and multiple structures.

The proposed method depends on two thresholds that control the flexibility of the model. Robust estimation of those parameters, in an adaptive scheme is an issue that need further elaboration and is a subject of ongoing research along with the application to other problems of methods [28, 45, 46] that could directly benefit from an effective solution to this problem.

Endnote

^a Matlab code available at http://www.cs.uoi.gr/~dgerogia