Sequential graph-based extraction of curvilinear structures

In this paper, a new approach is proposed to extract an ordered sequence of curvilinear structures in images, capturing the largest and most influential paths first and then progressively extracting smaller paths until a prespecified size is reached. The results are demonstrated both quantitatively and qualitatively using synthetic and real-world images. The method is shown to outperform comparator methods for certain cases of noise, object class, and scale, while remaining fundamentally easier to use due to its low parameter requirement.


Introduction
Extracting curvilinear structures, or sets of line segments, in digital images is an important low-level problem with many applications in computer vision [1].There have been many attempts to redefine the ideal centreline properties and to develop bespoke algorithms in a variety of domains, from medicine to engineering [1,2].
Existing methods introduce a wide range of parameters that are often difficult to tune and/or unintuitive.Furthermore, these parameters often require to be tuned by the user for each image, resulting in faults or robustness issues outside the chosen values [3].Other approaches [4] harness advances in machine learning to produce excellent results.The main problem with these approaches is that they require large training datasets, which are expensive in terms of data acquisition and model training, and may fail unpredictably in scenarios different from the training dataset.Training convolutional autoencoders and U-Net style architectures [5] gives excellent results in biomedical segmentation but generates blurry edges as a result of minimising a pixel-wise binary crossentropy loss function.
In this paper, we introduce a novel approach for sequential centreline extraction in two stages: (1) enhancement of the centreline based on prominent ridge detection and (2) extraction of the centreline based on an iterative graph-based optimisation.Graph-based centreline extraction approaches [6,7] are known to be robust to image noise, but they typically require user-specified start and end points.This approach searches through a subset of paths that fall along prominent ridges, sequentially returning the path with the largest intensity sum along its profile, and whose union with previous paths in the sequence maximally contributes to the output skeleton.Therefore, a user-specified seed is not required.The survey [8] indicates that simpler methods with fewer parameters lead to better generalisation without overfitting.The main contribution of the proposed method is that it only requires a single tunable parameter ε and two optimisation parameters t 1 and t 2 without initialisation regions or training data.The parameter ε corresponds to the length (in pixels) of the shortest lines to be captured.If the level of imaging noise is known a priori, ε is set to a slightly larger value than this level.The two optimisation parameters t 1 and t 2 are for the enhancement and extraction optimisation processes and do not require tuning on a per-image basis.Therefore, this proposed algorithm works remarkably well as an out-of-thebox solution.The main stages of our pipeline are shown in Fig. 1. 2 Related work

Centreline extraction
Many centreline extraction methods for 2D and 3D images have been proposed in the literature with applications in several fields [9].

Edge-based approaches
Simple gradient-based approaches are known to converge on local solutions, such as medial sheets in 3D.Gradient vector flow (GVF) addresses this problem by using a diffusion of the gradient information [2] in order to produce smooth centrelines.Leng et al. [10] propose enhancing the results of a simple Canny edge detector by using a voting field that enhances smooth curves and removes artificial lines caused by noise.Smistad [11] combined GVF, Hessian-based enhancement, and ridge-traversal.

Graph-based approaches
Graph-based approaches find the shortest path or separating cut between intensity or colour information.In particular, Rouchdy and Cohen [6] consider the accumulative overlap of shortest paths from the image boundary towards a user-specified seed, which is shown to be robust to image noise.Bauer [12] proposed a method to extract the centreline from airway images based on graph optimisation.Similarly, Jin et al. [7] find paths in objects with a tree-structure from a user-specified point, which is used to find the furthest point and the minimal connecting path.

Morphological-based approaches
These approaches are amongst the oldest in the field.[13] repeatedly applies morphological operations such as to shift boundary pixels without disturbing the object's local connectivity until it is a single pixel in width.Recently, [14] introduce a method to extract the curvilinear structure by applying the hit-or-miss transform in different directions followed by dilation operation to connect the broken lines.

Learning-based approaches
Recent learning-based extraction methods are more suitable to deal with the scene complexity problem in natural images [4,[15][16][17].Sironi et al. [4] method, learn huge numbers of parameters from large training datasets based on ground truths.As a result, their method is well suited to ill-posed problems, such as centreline extraction.Recent work by Shen et al. [15,17] proposed a multi-scaled learning framework that fuses the final output together.The output from their method is a binary image denoting the detected skeletons.

Surface-skeleton approaches
Traditional thinning approaches [18] iteratively remove outer layers of boundary pixels according to local stopping criteria.Similarly, geometric contraction approaches apply constrained iterative smoothing and/or merging, shrinking the surface into a thin centreline [19].However, these approaches are sensitive to surface noise.

Ridge detection
Ridges are defined as extrema of the image's largest surface curvature direction.Steger proposed a popular approach for extracting lines using Gaussian masks to estimate derivatives without bias in asymmetrical lines [3].Another approach by [20] handles crossings, junctions, and blobs by using anisotropic, multi-scale Gaussian kernels for second orderimage differentiation.

Overview
Enhancement is investigated based on the algorithmic definition of geographical prominence in 2D and 3D, which measures mountains from their summit to saddle points connecting to the next highest peak.Figure 2 shows how this property is desirable in images, whereby the most prominent peak (green arrow) captures the underlying object despite severe levels of multi-frequency inhomogeneous noise with higher local elevations, or summits, than the object itself.
The proposed method has two main conceptual stages: (1) centreline enhancement with geographical prominence and (2) sequential extraction of newly contributing bright paths.During the first step, an accumulative matrix is constructed using a voting procedure, whereby points in the matrix are at the index of prominent peaks of line (2D) or disc (3D) profiles.The second step searches through a subset of paths according to this matrix, sequentially extracting the most prominent paths that maximally contribute to the output skeleton according to a simple summation.

Enhancement process
In geography, the prominence of a peak in a mountain range is defined as the difference in elevation between the summit of the peak and the highest saddle that connects that summit to any higher area [21].This algorithm requires no parameters to be tuned (Table 1).
The proposed enhancement algorithm is formulated as a voting procedure based on prominence.Let us consider an image I (x) ∈ R n .If x is a random point in the image, n is a random unit vector, and l is a random length between 0 and the size of the smallest dimension of I , then a random intensity profile along a line from x to x + l n is defined as: (1) The intensity profile f is sampled and the index of the highest prominence peak is found, as in [21], as follows: and its pixel location: Then, a vote is cast into an accumulative map A: as shown in Fig. 3a.This process is repeated t 1 times, gradually accumulating more votes along ridge pixels.The accumulative map A is then normalised to be in the range [0, 1].

Extraction process
In order to extract the centreline, an undirected, weighted graph G(N , E) is constructed, where each pixel position x in the accumulative map A corresponds to a node in N .Graph edges E connect each node N at pixel x to its 8-connected neighbourhood of pixels.Graph node costs are defined by 1/A, such that the prominent ridges are set to small values and background regions tend to infinity.
From this graph definition G, a set of t 2 random shortest paths is produced iteratively: where each shortest path p i is calculated between two random nodes (according to the cost 1/A) using Dijkstra's algorithm [22].The most prominent path p is found and stored in P = { p} at each iteration, where p is defined: which is the segment of the path, from the set of paths in P, with the highest cost across the accumulative map A, that does not intersect paths found in previous iterations P. All paths in the set P must be disjointed (not sharing any pixels with other paths).Equation 6 is iterated (Fig. 4) until the length of the most prominent path p is shorter than ε.

Centreline in 3D images
The method is easily extensible to 3D images.In the enhancement procedure, the profile of a random disc is considered Fig. 4 The first three iterations of the sequential extraction process.After these three iterations, there is a large drop in the length of the next path, which is captured easily by ε instead of a line.In particular, a random plane is chosen, parametrised by a point x and normal n and the single most prominent peak over the set of lines (of length l) that sweep through 360 • on this plane is chosen.Only the most prominent peak from the set of lines on the plane will cast a vote in the accumulator array, see Fig. 3b.In the second step, the graph edges are defined as the full connectivity between the current pixel and its 26 neighbours instead of 8 neighbours Table 2 Comparison between the extraction centrelines obtained from the proposed method and the comparator methods for different noise types.The parameters for each method are tuned across all 6 scenarios and tuned manually per-image in the 2D version.The final skeleton is extracted as described in Sects.3.2 and 3.3.

Noise sensitivity
In particular, a 2D synthetic image of a curvilinear object with a known ground truth is created and the image is corrupted with different noise types, as shown in Table 2.The synthetic 2D curve is generated with a thickness radius of 3 and the resulting 2D image (128 × 128 pixels) is corrupted with different levels and types of noise.The Hausdorff distance (HD) is calculated between the extracted curve and its analytical ground truth (see Fig. 5).The results show that our method remains stable under severe quantities of noise, heavily corrupting the object to a PSNR of about 11.

Parameter space
In order to understand the effect of the optimisation parameters t 1 and t 2 in the two stages (Sects.3.2, 3.3), we conduct an Fig. 5 The HD between the extracted curve and the ground truth as shown in Table 2.These results use the same parameters: t 1 = 10000, t 2 = 2000, and ε = 2.The lines show the mean across 100 experiments and the error envelopes (transparent shaded regions) show the standard deviation experiment for 3 different synthetic images (128×128 pixels) and 3 different real-world images: road cracks (113 × 212), facial wrinkles (65 × 161), and fungal images (147 × 142).In particular, the mean Hausdorff distance (MHD) is just 3.5 pixels at the base of the hill (at the point t 1 = 8250, t 2 = 1950) for the synthetic images, and the MHD is Table 3 The HD measures corresponding to the images are shown in Table 2.The last two rows show the means and standard deviations of the HD across comparator methods

Comparison with the existing methods
The proposed method is also validated qualitatively against several 2D and 3D real-world images and the results are shown in Tables 5 and 6 , respectively.Images showing dark objects on bright backgrounds are inverted before processing.
In the 3D validation, the extracted curve is compared with the thinning method [26], which has a publicly available 3D implementation.
The proposed method is also validated against a wide range of centreline extraction methods using synthetic 2D data.In particular, we show (1) a modified traditional thinning approach [18], which has been improved with both median and Gaussian filtering and then thresholded to give a fair chance of success, (2) ridge detecting using Steger's popular approach of unbiased detector of curvilinear structures (UDCS) [3], (3) a Voronoi skeleton approach [25], and (4) the recent anisotropic multi-scale Gaussian kernels (AGK) [20]).
The results are shown in Table 2. Furthermore, the methods with their tuned parameters are compared across all 6 scenarios and are tuned manually on a per-image basis.In addition, the HD between our extracted curve and the ground truth is calculated for these 6 scenarios.The results are presented in Table 3.Our proposed approach succeeds under severe levels of noise where other approaches fail.The extracted centreline is robust, well-centred, and homotopic to the original object.However, it is not as smooth as with methods that use a Gaussian approximation of the derivatives.This can be addressed by applying 1D Gaussian smoothing in the Table 5 Results for different 2D images: (a) wires, (b) spirillum [27], (c) straight hair [28], (d) retina networks [29], (e-f) treponema [30,31], (g) telomeric DNA from HeLa cell clone [32], and (h) electron microscopic image of isolated mouse circular mtDNA [33] tangent direction for each line segment.In other approaches, the tuning of parameters requires significantly more effort.
In particular, the AGK method utilizes 11 different configurations to determine the set of parameter values.The UDCS method [3] gives excellent results when the parameters are carefully tuned on a per-scenario basis, but the method fails when the parameters are held constant.However, even simple approaches, such as thinning, outperform other methods when tuned on a per-image basis using different filtering techniques.In contrast, our method does not require such tuning.

Real-world validation
Our method was evaluated with the Ghent University Fungal Images dataset 1 (GUFI-1) [20], which is a popular and varied biological dataset with images of fungi growing in vitro.Each image has a resolution of 300 × 300 pixels and is accompanied with manually labelled ground truths for use in evaluating centreline extraction methods.These images contain a variety of ridges, different degrees of contamination, frequent overlaps and junctions, and large regions without edges and the dataset are therefore considered to be rather challenging.The ground truth centrelines are not always thin (with 1 pixel width), and therefore, the ground truth was manually altered to ensure the required thinness.The ground truth has been established independently by two experts, and the HD was then calculated between the two sets in order to ensure the validity of proposed ground truth.The MHD is 0.86 pixels, where the standard deviation is 0.54 pixels.The ground truths were revised by two other experts.The proposed method was evaluated alongside a variety of other methods, calculating the HD between the extracted skeletons and the ground truth.The results are shown in Fig. 6 and Table 4.The results show that our method comes closest to the ground truth across the dataset; however, certain features are missing which may be better captured by other methods [3,18,20].However, unlike other methods, our approach tends to be cleaner and results in fewer artifacts (i.e.spurious branches).It is inferred that the proposed method is beneficial in the task of collecting reliable metrics across a large dataset, but, for some specific analyses of individual images, it is recommended to consider more sensitive methods, such as [3] (Tables 5, 6).
Table 6 Comparison between our method and topological thinning [26] using 3D images of: neocortical axon trees (a) [34], olfactory projection neuron tree (b) [34], and visual cortical trees (c) [34].The left columns show top-down views with the maximum intensity projection of the 3D object structure

Future work
In our implementation,1 the performance of the algorithm may be significantly improved by removing low-weighted edges in the graph according to a thresholding parameter with very little impact on quality.Indeed, this is especially useful for segmenting large 3D images.In this paper, we have explored the performance of the proposed algorithm and the optimisation of our method is left as an area of future work.
Although the proposed method gives good 2D/3D results for biomedical and non-biomedical images, there is room for improvement.In particular, the enhancement process should be removed and the graph cost function should be extended in order to ensure centeredness in flat intensity distributions, as well as the better handling of junctions.This extension could also consider a distance metric that gives preference

Conclusion
A novel approach is described which searches through a subset of paths that fall along prominent ridges, sequentially returning the most prominent path that maximally contributes to the output centreline.The method is evaluated against real-world images and is shown to be comparable to the comparator methods with extracted centrelines that are close to the ground truth.In particular, the method has different strengths and weaknesses, such as its ability to bridge inhomogeneous gaps (Table 6b) and to handle complex/busy backgrounds.Furthermore, the method requires little tuning and works well as an out-of-the-box solution in medical (e.g.blood-vessel extraction) and biological (e.g.plant roots and neural networks) imaging applications.

Fig. 1
Fig.1Workflow of the proposed approach: a t 1 intensity profiles, defined by random lines parametrised with random positions x directions n and lengths l, are considered and cast votes at the locations of the maximum peaks into an accumulative space (b).This is followed by defining a graph G(N , E) used to build a set of t 2 shortest paths

Fig. 2
Fig. 2 An intensity profile along a random line.a The area of interest in image has inhomogeneous backgrounds, highlighted by red dashed square and the random line in blue, and b profile that shows prominent peaks compared with the most prominent peak (green arrow) that captures the underlying object structure (colour figure online)

Table 1 Fig. 3
Fig. 3 An illustration of finding the most prominent peaks in 2D and 3D images.a The proposed method uses a random line in 2D to find the most prominent peak x and accumulates it into A. b The most prominent peak on a random plane in 3D

Fig. 6 A
Fig. 6 A selection of images in the Ghent University Fungal Images (GUFI-1) dataset [20] alongside the ground truth centrelines and extracted curves from multiple methods.a Original image, b ground truth, c thinning [18], d UDCS [3], e Voronoi [25], f AGK [20], and g our method

Table 4
Comparing the proposed method with the comparator methods in terms of the HD using the GUFI-1 dataset