1 Introduction
Many computer vision algorithms have been proposed to analyze treelike curvilinear structures (also called line networks), such as blood vessels for retinal images
[7, 27], filament structure of biological images [24, 33], road network for remote sensing [11, 19, 32], and defects on materials [5, 12]. Although human can intuitively perceive the curvilinear structures in very different images and application scenarios, most curvilinear structure detection methods work only for one specific application. The difficulties arise because the geometry of the line networks shows various shapes. Also, there is often insignificant contrast between a pixel belonging to a curvilinear structure and the background texture. Thus, the information from an individual pixel fails to interpret the such topology.methods, (d) the proposed algorithm represents well topological features of the curvilinear structure. In this example, road network is partially occluded by trees or cars. Setting a threshold value for the local measure often fails to detect the underlying curvilinear structure because this process neglects correlated information of the reconstructed curvilinear structure. Although the centerline is able to quantify scale (width) of curvilinear structure, it inaccurately classifies pixels around intersections. In this work we learn structured output ranking of the image patches that contain a part of curvilinear structures. The proposed graphbased representation uses the smallest pixels but also preserves the topological information.
To represent an arbitrary shape of the curvilinear structure, segmentation algorithms [2, 7, 20] evaluate linearity scores, which correspond to likelihood of pixels given the curvilinear structure. Then, a threshold value is set to discard pixels showing low scores. This process implicitly ignores correlated information of the pixels on the curvilinear structure, so that it yields discontinuities on the reconstruction results. Centerline detection algorithm [26] is able to encode the width of the curvilinear structure; however, it is weak for localizing junctions of the curvilinear structure. On the other hand, graphbased models [9, 24, 33, 34] reconstruct the curvilinear structure with geodesic paths in a sparse graph. For an efficient computation, these algorithms usually exploit a subset of pixels which correspond to local maxima of the linearity score. Contrary to the previous methods, our approach initially infers the curvilinear structure based on the structured ranking scores, which are related to the spatial patterns of the latent curvilinear structure. We then search for the longest geodesic path on the graph and iteratively add fine branches to rebuild the curvilinear structure with minimal pixels. Figure 1 compares the results obtained with curvilinear structure segmentation [2], centerline detection [26], and the proposed method.
In this work we aim to learn structured rankings of the input image patches which may contain a part of curvilinear structures. In particular, we assign a high ranking score for the image patch if it shows plausible spatial pattern to the latent curvilinear structures. Since the curvilinear structures are arbitrarily oriented within in an image patch, we measure the local orientation of the image patch and rotate it with respect to the baseline orientation (align horizontally). We then explore geodesic paths in the graph built upon the structured ranking score map. The proposed algorithm can provide the topological importance level of the curvilinear structure. Our experiments demonstrate that the proposed algorithm localizes the curvilinear structure with a high accuracy compared to the stateoftheart methods.
1.1 Related work
Geometrically speaking, curvilinear structures are rotatable and elongated shape. Image gradient information [7, 8, 13, 20, 25]
is useful to estimate local orientation and to describe the shape of the curvilinear structure. The local image features are insufficient to discern curvilinear structures from undesirable highfrequency components due to the lack of geometric interpretation. To take the structural information into account, graphical models
[9, 30] formulate an energy optimization problem with spatial constraints in a local configuration. Stochastic models [14, 19] also specify a distribution of the line objects given image data with pairwise interaction terms. A sampling technique [10]is employed to maximize the probability density. However, the geometric constraints are heuristically designed, and the number of constraints are increased to describe complex shaped line networks. Recently, machine learning algorithms have been involved to detect curvilinear structures latent in various types of images. Becker
[2] applied a boosting algorithm to obtain an optimal set of convolution filter banks. A regression model is proposed to detect centerlines by learning the scale (width) of the tubular structures with the nonmaximum suppression technique [26].Structured learning system has been employed in image segmentation models based on random fields [3, 16, 21, 28]. More specifically, Structured Support Vector Machine (SSVM) [29] is used to predict model parameters for inference of the structured information between input image space and output label space. Exploring all possible combinations of the labels in the output space is computationally intractable. Thus, the random field models define the pairwise relationship in a neighborhood system to enforce the labeling consistency. While such prior models based on random fields are successful to describe convex shaped objects, it is inefficient to detect thin and elongated shaped curvilinear objects.
The main contributions of this work are summarized as follows:

We propose an orientationaware curvilinear feature descriptor for the curvilinear structure inference;

We learn a ranking function which can infer spatial patterns of the curvilinear structure within image patches; and

We reconstruct the latent curvilinear structure based on graph topology with the structured output rankings.
The rest of the paper is organized as follows: Sec. 2 provides the mathematical description and outline of our method. Sec. 3 proposes an orientationaware curvilinear feature descriptor based on oriented image gradients. Sec. 4 explains how we infer the structured output rankings of the given input image patches. Sec. 5 develops a graphbased curvilinear structure reconstruction algorithm. Sec. 6 shows experimental results on different types of datasets. Finally, Sec. 7 concludes this work.
B 
2 Overview
In this section, we define notations and provide an overview of the proposed algorithm (see Figure 2). Assume that an image contains a curvilinear structure. We denote the latent curvilinear structure for any pixel of the image :
(1) 
This function is also related to a ground truth map, which is manually labeled, for the machine learning framework and for the performance evaluation. We compute a curvilinear feature map that represents oriented gradient information (see Section 3.1). Since information embedded in a single pixel is limited to infer the latent spatial patterns, we exploit image patches to compute input feature vectors. Let be a patch of the feature map values within size of square window centered at , , . Using a rotation matrix defined on Euclidean image space, we can rotate a patch with respect to the given orientation such as .
We use graph theory for shape simplification of the complex curvilinear structure. For the pixels showing higher rankings, we build an undirected and weighted graph . Then, we look for the longest geodesic path which corresponds to the coarset curvilinear structure in the image. We iteratively reconstruct the curvilinear by collecting paths that connect the remotest vertices in the graph. Consequently, the proposed algorithm can represent the different levels of detail in the latent curvilinear structure using the minimum number of pixels.
3 Curvilinear Feature
This section is devoted to compute the curvilinear feature descriptor that is used for the inputs of the learning system. We perceive the latent curvilinear structure based on inconsistency of background textures and its geometric characteristics. In other words, a sequence of pixels corresponding to the curvilinear structure has different intensity values compared to its surroundings, and shows thin and elongated shape. Thus, we compute multidirection and multiscale image gradients to detect locally oriented image features.
3.1 Curvilinear feature extraction
We obtain oriented gradient maps using a set of steerable filters [8, 13, 25]. Before applying the convolution operations, we normalize the training and test images to remove the effects of various illumination factors:
(2) 
where is the sample mean of the image.
The steerable filters are created by the second order derivative of isotropic 2D Gaussian kernels. Let be a steerable filter that accentuates image gradient magnitude for direction at scale . To take into account varying orientations and widths of the curvilinear structure, a feature map is combined by multiple filtering responses:
(3) 
where denote the total number of orientations. We first find the maximum filtering responses of the scale spaces, and then average for all directional responses. In this work we consider 8 orientations and 3 scales.
3.2 Local orientation estimation
In the previous subsection, we evaluate the presence of curvilinear structure from the amount of gradient magnitudes in image patches. Complex shaped curvilinear structure also consists of various local orientations.
Assume that be a basis spatial pattern of the curvilinear structure for the baseline orientation (). We manually design B as a simple line shape template which highlights the middle of image patch (see Figure 3). If an image patch contains a part of curvilinear structure, pixels on the curvilinear structure are sequentially accentuated towards a particular direction . Recall that the geometric properties of the curvilinear structure, it is rotatable and symmetric. We can rotate image patch to be aligned with the basis spatial pattern.
More specfically, we estimate the local orientation shown in image patches using statistical difference of the areas determined by B. For a rotated image patch , we compute normalized histograms of the pixels which are labeled as curvilinear structure and its counterpart , where denotes histgram of the given distribution with bins. Similar to [1], we employ test to compute statistical difference of two distributions:
(4)  
(5) 
Figure 3 compares statistics of the image patches with respect to different orientations. If an image patch contains curvilinear structure, test score shows unimodal distribution which shows a peak at the dominant orientation.
4 Learning
We regulate the orientation and the shape of the input patches with respect to the basis spatial pattern B. In this section, we aim to learn a function that predicts structured output rankings of the input image patches. Thus, the output rankings indirectly infer the spatial patterns of the curvilinear structure for image patch.
Let be an operator to convert an image patch into a column vector. For a pixel , the input feature vector is defined as and denote the corresponding rank value. For a feature vector, we exploit subsampled pixels on the linear template to avoid data imbalance. Thus, the dimension of feature vector is small than the total number of pixels within the patch, , .
For the setup of machine learning, a training dataset consists of the inputandoutput pairs. Let be an unordered list of the input feature vectors and be the corresponding classes. The input feature vector z is built from the curvilinear feature map . Our goal is to learn a ranking function which assigns a global ordering (ranking) of feature vectors: . The inner product of the model parameter and a feature vector is used to predict ranking score of the given image patch.
In our work Structured SVM framework [15, 23] is employed to exploit the structure and dependencies within the output space . To encode the structure of the output space, a loss value of each input feature vector is defined: . The loss value evaluates the quality of the learning system. Intuitively, the input image patches containing curvilinear structure are comparable to the basis spatial pattern. Hence, as an image patch is similar to the basis spatial pattern, we should put it on a higher ranking. Specifically, we compute the loss value as the overlapping ratio between image patch from the groundtruth map and the basis :
(6) 
where denotes a rotated image patch referring to groundtruth map . The objective function of Structured SVM with a single slack variable is given by:
(7) 
where is the indicator variable to reduce the number constraint in linear complexity :
(8) 
Cutting plain algorithm [29] is employed to solve (7). For the details about the Structured SVM optimization, we refer the reader to [15, 23, 29]. In the following section, we plot the initial segmentation map from the structured rankings scores. We also exploit this ranking scores for the shape simplification based on graph traversal algorithm.
5 Curvilinear Structure Reconstruction
To represent the latent curvilinear structure, various local measures have been proposed to detect irregular shaped curvilinear structure. The score of the local measure is regarded as a likelihood probability whether a pixel is on the latent curvilinear structure. Most of the previous works represent the curvilinear structure as a binary map based on the linearity scores of the pixels. This approach easily misinterpret the topological features. In this section, we propose a novel structured score map based on the output ranking scores and the basis spatial pattern. We also develop a graphbased model which is able to organize the topological features of the structure in different levels of detail. Figure 4 shows stepbystep processing results to reconstruct the latent curvilinear structure.
5.1 Structured score map
We have obtained a score function that evaluates the compatibility between the input features and the underlying curvilinear structure based on the structured output rankings. From the training dataset, we precompute the proportion of pixels being part of the curvilinear structure. Note that the value maximizes test scores at the groundtruth maps in the training dataset. During the test phase, we retain the output ranking according to the output ranking scores from the highest to the th rank. As shown in Figure 4 (c), a binary map can be disconnected.
To avoid unwilling breakpoints, we composite the structured score map using the binary map, induced by ranking scores , and the basis spatial pattern B. Recall that we initially estimate the ranking score using the shape similarity between oriented image patch of ground truth and the basis spatial pattern in (6). Thus, it can be regarded as the inverse mapping from output ranking score to the latent curvilinear structure. Let be the column vector version of the binary mask B and be the patch of structured score map centered at with . We minimize the following cost function with respect to :
(9) 
This is a least square problem, so that the solution is found at points which satisfies . To reconstruct the structured segmentation map, we synthesize obtained patches at the subsampled grid points on image , where . We refer the readers to [18] for the implementation details of the related texture synthesis technique.
5.2 Progressive curvilinear path reconstruction
We consider the graph where is the set of pixels. Two pixels and are connected by the edge if and only if . Moreover, we assign a weight for each edge as . A path in the graph is a sequence of distinct vertices such that consecutive vertices are adjacent. The length of is the sum of the weights of its edges, and the distance between two vertices and is the minimum length of a path from to . The eccentricity denotes the maximum distance from the vertex to a vertex , , . The diameter of equals , , it is the maximum distance between two vertices in .
To simplify the latent curvilinear structure, we look for long geodesic paths in the subgraph of induced by the pixels with structured score map . More precisely, our algorithm computes a diameter of , , a shortest path with length . This path is added in the simplified curvilinear structure, then the path is contracted into a single (virtual) vertex. For the implementation, the weight of all edges of become . We repeat this process till the diameter of the subgraph is larger than predefined path length . The entire procedure is summarized in Algorithm 1.
The timeconsuming part of the proposed algorithm is the computation of a diameter of . Rather than computing all pair distances (which requires a linear number of application of Dijkstra’s algorithm), we use an efficient heuristic algorithm called 2sweep algorithm [6]. The 2sweep algorithm randomly picks a vertex in , then performs Dijkstra’s algorithm from to find a vertex at the maximum distance from , , . Then, it computes (using Dijkstra’s algorithm) a path from to a vertex at the maximum distance from . The length of the second path (from to ) is a good estimation of . Note that this algorithm is able to compute the exact diameter if the graph has tree structure [4]. Topologically speaking, most of latent curvilinear structures are very close to trees [9, 31]. Therefore, the 2sweep algorithm is well adapted to reconstruct the treelike curvilinear structures. For a better understanding, we schematically explain the intermediate steps of the proposed curvilinear structure simplification algorithm in Figure 5.
6 Experimental results
In this section, we first discuss the parameters of the proposed algorithm and datasets. We then compare the quantitative and qualitative results of the proposed algorithm and those of competing models proposed by [7], [20], [2], and [26].
6.1 Parameters and Datasets
The proposed algorithm requires few parameters to compute the curvilinear feature descriptor . We use 8different orientations in this work: . The size of steerable filters is fixed to pixels. The scale factors and the minimum path length are adaptively selected for each dataset to obtain the best performances. The binary spatial pattern B is based on the physical attribute of the dataset in that thickness of curvilinear structure shows various depending on the dataset. For Structured SVM training, we sample 2000 image patches on each dataset. All patch sizes are fixed to . which controls the relative importance of slack variables is set to for all datasets. We choose the set of parameters for the SSVM training via 3fold cross validation [17] which maximizes the average score [22] of the training set.
We test our curvilinear structure model on the following public datasets:

DRIVE [27]: The dataset consists of 40 retina scan images with manual segmentation by ophthalmologists to evaluate the blood vessel segmentation algorithms. We use 20 images for the training and 20 images for the test, respectively. The path length is set to . The scale factors and thickness are set to and , respectively.

RecA [14]: We collect electron microscopic images of RecA proteins on DNA which contain filament structure. We use 4 training images and 4 test images. We use , , and .

Aerial [26]: The dataset contains 14 remote sensing images of road networks. We select 7 images for the training and 7 images for the test, respectively. We use , , and .

Cracks [5]: Images of the dataset correspond to road cracks on the asphalt surfaces. We use 6 images to train and test the algorithms on different 6 images. For this dataset, we set to , , and , respectively.
6.2 Evaluations
The proposed algorithm progressively reconstructs the curvilinear structure by adding a long path on the subgraph. Figure 6 shows the intermediate steps of the proposed curvilinear structure simplification algorithm for DRIVE dataset. Unlike the previous models, the proposed algorithm is able to show different levels of detail for the latent curvilinear structure. Such information to visualize shape complexity of the curvilinear structure cannot be retrieved by setting a threshold. In practice, a few number of iterations is required to converge the algorithm and each step to find a long path takes less than milliseconds for the computation. For the experiments, we use a PC with a 2.9 GHz CPU (4 cores) and 8 GB RAM. Moreover, we visually compare the performance of the proposed algorithm with the competing algorithms. Figure 7 shows the results from DRIVE, RecA, Aerial, and Cracks dataset, respectively. The proposed algorithm is the most suitable to show the topological information of the latent curvilinear structures.
For the quantitative evaluation, we provide scores of the proposed algorithm and the stateoftheart models in Table. 1. The measure of true positive is sensitive for the misalignment; therefore, we consider surrounding pixels of the detection results as the true positive if a predicted point is falling into the ground truth with a small radius (equivalent to its thickness parameter ) similarly to [26]. We also provide the average proportion of pixels to represent the curvilinear structures. It shows that the proposed algorithm efficiently draws the curvilinear structures using smaller number of pixels than the other algorithms.
The proposed algorithm achieved the best scores for all datasets except the DIRVE dataset. It is because the proposed algorithm use the fixed thickness to describe linear structure in the binary mask ; however, the images consisting of DRIVE dataset exhibit varied thicknesses of blood vessels and many junction points. In other words, we obtain a model parameter regarding to the manually designed binary pattern . To overcome this drawback, we plan to exploit generative binary patterns based on the training images in the future, which is possible in our framework.
DRIVE  RecA  Aerial  Cracks  

Frangi [7]  0.33  0.33  0.32  0.056  
Law & Chung [20]  0.43  0.21  0.25  0.085  
(a)  Becker [2]  0.50  0.45  0.53  0.23 
Sironi [26]  0.55  0.50  0.55  0.27  
Proposed (Graph)  0.36  0.59  0.59  0.38  
Frangi [7]  37.09  16.14  22.19  28.64  
Law & Chung [20]  19.60  29.82  29.51  33.00  
(b)  Becker [2]  12.67  9.39  9.76  32.13 
Sironi [26]  5.41  5.47  0.83  17.12  
Proposed (Graph)  0.17  0.57  0.35  1.07 
The topology of the latent curvilinear structure is varied while we compare it within the same dataset. Especially, Cracks dataset consists of difficult images due to the rough surfaced background textures. The normalization operation (2) is necessary to remove irregular illumination factors on the datasets. Manual segmentation (ground truth) contains many errors around boundaries and minutiae components. Also, the number of training data employed in this work is relatively small due to the difficulties of making the accurate annotations. It is remarkable that the proposed algorithm shows good performance for all datasets using the imperfect training dataset.
7 Conclusions
This paper proposed a curvilinear structure reconstruction algorithm based on the ranking learning system and graph theory. The output rankings of the image patches corresponded to the plausibility of the latent curvilinear structure. Using an optimal number of pixels, the proposed algorithm provided different levels of detail during reconstruction of the curvilinear structure. More precisely, we learned a ranking function based on Structured SVM with the proposed orientationaware curvilinear feature descriptor. In this paper, we developed a novel graphical model that infers the curvilinear structure according to the topological importance. The proposed algorithm looked for remote vertices on the subgraph which is induced from the output rankings. Across the various types of datasets, our model showed good performances to reconstruct the latent curvilinear structure with a smaller number of pixels comparing to the stateoftheart algorithms.
References
 [1] P. Arbeáez, M. Maire, C. Fowlkes, and J. Malik. Contour detection and hierarchical image segmentation. IEEE TPAMI, 33(5):898–916, 2011.
 [2] C. Becker, R. Rigamonti, V. Lepetit, and P. Fua. Supervised feature learning for curvilinear structure segmentation. In MICCAI, pages 526–533, 2013.
 [3] L. Bertelli, T. Yu, D. Vu, and B. Gokturk. Kernelized structural SVM learning for supervised object segmentation. In CVPR, pages 2153–2160, 2011.
 [4] M. Borassi, P. Crescenzi, M. Habib, W. A. Kosters, A. Marino, and F. W. Takes. Fast diameter and radius BFSbased computation in (weakly connected) realworld graphs: With an application to the six degrees of separation games. Theor. Comput. Sci., 586:59–80, 2015.

[5]
S. Chambon, C. Gourraud, J.M. Moliard, and P. Nicolle.
Road crack extraction with adapted filtering and Markov modelbased segmentation.
In VISAPP(2), pages 81–90, 2010.  [6] D. G. Corneil, F. F. Dragan, M. Habib, and C. Paul. Diameter determination on restricted graph families. Discrete Applied Mathematics, 113(2–3):143–166, 2001.
 [7] A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever. Multiscale vessel enhancement filtering. In MICCAI, pages 130–137, 1998.
 [8] W. T. Freeman and E. H. Adelson. The design and use of steerable filters. IEEE TPAMI, 13(9):891–906, 1991.
 [9] G. González, E. Türetken, F. Fleuret, and P. Fua. Delineating trees in noisy 2D images and 3D imagestacks. In CVPR, pages 2799–2806, 2010.

[10]
P. J. Green.
Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.
Biometrika, 82(4):771–732, 1995.  [11] J. Hu, A. Razdan, J. C. Femiani, M. Cui, and P. Wonka. Road network extraction and intersection detection from aerial images by tracking road footprints. IEEE TGRS, 45(12):4144–4157, 2007.
 [12] S. Iyer and S. Sinha. A robust approach for automatic detection and segmentation of cracks in underground pipeline images. Image and Vision Computing, 23(10):921–933, 2005.
 [13] M. Jacob and M. Unser. Design of steerable filters for feature detection using Cannylike criteria. IEEE TPAMI, 26(8):1007–1019, 2004.
 [14] S.G. Jeong, Y. Tarabalka, and J. Zerubia. Marked point process model for curvilinear structures extraction. In EMMCVPR 2015, LNCS 8932, pages 436–449, 2015.
 [15] T. Joachims. Training linear SVMs in linear time. In KDD, pages 217–226, 2006.
 [16] S. Kim, C. D. Yoo, S. Nowozin, and P. Kohli. Image segmentation using higherorder correlation clustering. IEEE TPAMI, 36(9):1761–1774, 2014.
 [17] R. Kohavi. A study of crossvalidation and bootstrap for accuracy estimation and model selection. In IJCAI, pages 1137–1143, 1995.
 [18] V. Kwatra, I. Essa, A. Bobick, and N. Kwatra. Texture optimization for examplebased synthesis. In SIGGRAPH, pages 795–802, 2005.
 [19] C. Lacoste, X. Descombes, and J. Zerubia. Point processes for unsupervised line network extraction in remote sensing. IEEE TPAMI, 27(10):1568–1579, 2005.
 [20] M. W. Law and A. Chung. Three dimensional curvilinear structure detection using optimally oriented flux. In ECCV, pages 368–382, 2008.
 [21] A. Lucchi, Y. Li, and P. Fua. Learning for structured prediction using approximate subgradient descent with working sets. In CVPR, pages 1987–1994, 2013.
 [22] D. R. Martin, C. C. Fowlkes, and J. Malik. Learning to detect natural image boundaries using local brightness, color, and texture cues. IEEE TPAMI, 26(5):530–549, 2004.
 [23] A. Mittal, M. B. Blaschko, A. Zisserman, and P. H. S. Torr. Taxonomic multiclass prediction and person layout using efficient structured ranking. In ECCV, pages 245–258, 2012.

[24]
H. Peng, F. Long, and G. Myers.
Automatic 3D neuron tracing using allpath pruning.
Bioinformatics, 27(13):239–247, 2011.  [25] P. Perona. Deformable kernels for early vision. IEEE TPAMI, 17(5):488–499, 1995.
 [26] A. Sironi, V. Lepetit, and P. Fua. Multiscale centerline detection by learning a scalespace distance transform. In CVPR, pages 2697–2704, 2014.
 [27] J. J. Staal, M. D. Abramoff, M. Niemeijer, M. A. Viergever, and B. van Ginneken. Ridge based vessel segmentation in color images of the retina. IEEE TMI, 23(4):501–509, 2004.
 [28] M. Szummer, P. Kohli, and D. Hoiem. Learning CRF using graph cuts. In ECCV, pages 582–592, 2008.
 [29] I. Tsochantaridis, T. Joachims, T. Hofmann, and Y. Altun. Large margin methods for structured and interdependent output variables. Journal of Machine Learning Research (JMLR), 6:1453–1484, 2005.
 [30] E. Türetken, F. Benmansour, B. Andres, H. Pfister, and P. Fua. Reconstructing loopy curvilinear structures using integer programming. In CVPR, pages 1822–1829, 2013.
 [31] E. Türetken, G. González, C. Blum, and P. Fua. Automated reconstruction of dendritic and axonal tress by global optimization with geometric priors. Neuroinformatics, 9(2–3):279–302, 2011.
 [32] S. Valero, J. Chanussot, J. A. Bendiktsson, H. Talbot, and B. Waske. Advanced directional mathematical morphology for the detection of the road network in very high resolution remote sensing images. Pattern Recognition Lett., 31(10):1120–1127, 2010.
 [33] Y. Wang, A. Narayanaswamy, and B. Roysam. Novel 4D opencurve active contour and curve completion approach for automated tree structure extraction. In CVPR, pages 1105–1112, 2011.
 [34] T. Zhao, J. Xie, F. Amat, N. Clack, P. Ahammad, H. Peng, F. Long, and E. Myers. Automated reconstruction of neuronal morphology based on local geometrical and global structural models. Neuroinformatics, 9(2–3):247–261, 2011.
Comments
There are no comments yet.