Curvilinear Structure Modeling: Reviews

Keywords: curvilinear structure modeling, marked point process, structure inference, machine learning Résumé modélisation de structure curviligne, processus ponctuel marqué, inférence de structure, apprentissage automatique reconstruction . . . .

In this dissertation, we propose curvilinear structure reconstruction models based on stochastic modeling and ranking learning system. We assume that the entire line network can be decomposed into a set of line segments with variable lengths and orientations. This assumption enables us to reconstruct arbitrary shapes of curvilinear structure for different types of datasets. We compute curvilinear feature descriptors based on the image gradient profiles and the morphological profiles. For the stochastic model, we propose prior constraints that define the spatial interaction of line segments. To obtain an optimal configuration corresponding to the latent curvilinear structure, we combine multiple line hypotheses which are computed by MCMC sampling with different parameter sets. Moreover, we learn a ranking function which predicts the correspondence of the given line segment and the latent curvilinear structures. A novel graph-based method is proposed to infer the underlying curvilinear structure using the output rankings of the line segments. We apply our models to analyze curvilinear structure on static images. Experimental results on wide types of datasets demonstrate that the proposed curvilinear structure modeling outperforms the state-of-the-art techniques.

Acronyms and Notations

Acronyms

Image gradient vector at pixel x i (σ 2 is for computing derivatives) H(x i ; σ 2 ) Hessian matrix at pixel x i (σ 2 is for computing derivatives)

Stochastic Process

K(X, X ) Transition kernel which updates state from X to X K m (X, X ) Sub-transition kernel to support m-type jump from X to X q m Probability to select

K m (•, •) Q m (X |X)
Proposal probability to support m-type jump from X to X α m (X, X )

Acceptance ratio to support m-type jump from X to X

Graph Theory

G Undirected graph V Vertex set, V = {v} E Edge set, E = {e ij = (v i , v j )} |V |
The number of vertices in V dist (u, v) The minimum length of path from u to v in the graph G ecc (v) Eccentricity of vertex v (i.e., ecc(v) = max u∈V dist(u, v)) diam (G) Diameter of graph G (i.e., diam(G) = max v∈V ecc(v))

Machine Learning

2.3

We compare convolution filter banks which are built constructively (top row) and analytically (bottom row). Increasing smoothing factor σ is related to estimating width of the structure, it is more likely to remove noises. On the other hand, controlling the aperture value a of constructive filters provides consistent result that highlight pixels on the curvilinear structure according to the width. The associated convolution filters are depicted in left bottom of each image. . . . . . . . . . . . .

2.4

Toy example to show how the path opening operator finds valid paths on the adjacency graph: (a) an adcaency graph composed with south-north cones, (b) greed nodes are associated with pixels after image binarization, (c) seeking paths, denoted red color, on the graph, and (d) a path is valid if its length is larger than pre-defined length L. . . . . . . . . .

Morphological profiles:

Path opening [ Talbot 07] results with different length parameter L. As the length parameter is increased, the path opening operator removes small objects which are not on a valid path. 3.6 To prevent congestion of line segments, we compute the overlapping area of pairwise line segments. We take into account either pixels on the line segment s m and its surroundings A(s m ), which are denoted by gray color in this figure. The criteria function Υ(s m , s m ) computes the number of pixels falling into the overlapping area of two line segments. We allow small overlapping between two line segments (a), whereas we reject line configuration (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.7

The latent curvilinear structure, denoted by gray color, is reconstructed by a set of line segments. A state of the Markov chain is related to a line configuration, and updates its state according to the transition kernels. We finally obtain an equilibrium state of the Markov chain which maximizes the corresponding probability density in the sampling space. . Mathematical morphology operator, path opening [ Talbot 07], is applied on such gradient magnitude image (c). Linearity score of each pixel is drawn by the supervised feature learning algorithm [Becker 13] (d). We provide line hypotheses (e)-(g) associated with different hyperparameter vectors. Composition result (h) is equivalent to mixture probability density, and it highlights pixels corresponding to linear structures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.11

We provide intermediate sampling processes when the temperature parameter T t is decreasing. The results shown in first row are obtained without specifying seed segment. For the second row, we randomly set 20 seed segments and run the algorithm. For the third row, we initialize 20 line segments which are highly corresponding to underlying curvilinear structures. The algorithm converges toward almost the same solution regardless of the initial state. . . . . . . . . . . . . . . . . . .

3.

12 Precision-and-recall curves for pixelwise segmentation of curvilinear structures using path opening operator [ Talbot 07] with different setups of length, supervised feature learning [Becker 13], baseline MPP, and the proposed method. . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.13

We visualize the localization of the curvilinear structures on a Retina image. Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to ground truth (b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground truth. . . . . . . . . . . . . . . . . . . . . .

3.14

We visualize the localization of the curvilinear structures on a DNA image. Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to ground turth(b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground truth. . . . . . . . . . . . . . . . . . . . . .

3.15

We visualize the localization of the curvilinear structures on a wrinkle image. Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to ground truth (b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground-truth. . . . . . . . . . . . . . . . . . . . . .

3.16

We visualize the localization of the curvilinear structures on a Crack image. Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to ground truth(b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground truth. . . . . . . . . . . . . . . . . . . . . .

4.1

Comparing with the segmentation ((c), [Becker 13]) and the centerline detection ((d), [Sironi 14]) methods, the proposed algorithm is able to organize the topological importance of the curvilinear structure. The output rankings (e) highlight areas according to the geometric similarity of line segments and the underlying curvilinear structure. Our graphical model represents different levels of detail of the latent curvilinear structure with the minimum number of pixels. The proposed algorithm reconstructs the coarsest structure (f), and then iteratively adds fine branches (g). The final result (h) is obtained by a few iterations. . . .

4.2

For the same data points, we compare two linear separators. Intuitively, the separator in the bottom row is safer in that its decision boundary can tolerate a large noise. The maximum noise tolerance of a separator is called its margin corresponding to the thickness of the gray areas in this figure. Thus, SVM is designed to find a hypothesis with the maximum margin for the given training data points. . . . . . . . . . . . . . . . .

4.4

For the given linearly non-separable data points, we employ soft margin SVM which allows small classification errors around margin. A slack variable measures the such violation on the marginal areas. We modify the objective function of SVM to take into account the amount of margin violation with a parameter C. The hyperplane on the left figure is obtained using a small C so that it is less sensitive for classification errors and margin violations. On the other hand, the hyperplane on the right figure is more strict with respect to the margin violation by employing a larger C value. . .

Introduction

Introduction

F IELD of research and industry related to computer vision has been developed rapidly recent years due to the growth of computational power and image capturing devices. Now various types of image data are available for analysis. From microscopic images to satellite data, we observe curvilinear structures as different forms. It is required to detect the latent curvilinear structures to improve visibility, to infer biological information, and to extract geographical information. While human can intuitively distinguish these latent curvilinear structures, a technical definition of the curvilinear structure is unclear to date. The difficulties arise because the topology of the curvilinear structures exhibits varied shapes; thus, it is impossible to describe complex curvilinear structure with a few parameters. Also, depending on the acquisition environment, the curvilinear structures can be hidden under the background textures with a low contrast. Therefore, many data representation and computer vision algorithms to detect curvilinear structures have been developed for specific applications. Figure 1.1 shows some curvilinear structures for different applications. In early vision, researchers [Freeman 91, Perona 95, Jacob 04] have developed image convolution filters to decompose curvilinear features from background textures. The main idea behind such filter design is to create line shape templates to extract locally oriented gradient information. Frangi et al. analyzes the eigenvalues of the Hessian matrix for a given image to obtain principle directions of the local structure [ Frangi 98]. Optimally Oriented Flux (OOF) measures the amount of outgoing gradient flux to find the curvilinear structures [ Law 08]. Morphological operator collects pixels according to the structural similarity on the elongated path [ Talbot 07]. However, due to the lack of shape interpretation, such methods based on local image features are insufficient to reconstruct underlying curvilinear structure. On the other hand, graphical models such as [González 10,Türetken 13b] define geometric constraints in a local configuration and globally minimize their cost function. More precisely, the graph-based algorithms initialize some points highly corresponding to the latent curvilinear structure, and then define a path which connects these points with geometric priors to provide plausible shapes. Geometric properties of line network are involved as constraint terms when an energy optimization problem is formulated. Stochastic models [Lacoste 05, Jeong 15a] reconstruct curvilinear structures by sampling multiple line segments to maximize a posterior probability of given image data. Similarly to the graph-based representation, geometric priors are considered to define the connectivity and curvature of line segments. Recently, machine learning algorithms have been proposed to detect curvilinear structure. Becker et al. [Becker 13] applied a boosting algorithm to obtain an optimal set of convolution filter banks. Sironi et al. [Sironi 14] developed a regression model to estimate the scale (width) of the curvilinear structures and to localize the centerlines.

Although the contour grouping algorithms [Tu 06, Arbeáez 11] examine image features corresponding to curves and lines, the goal is quite different from the curvilinear structure reconstruction techniques. The contour grouping algorithms seek closed contour lines to divide an image into meaningful regions. Therefore, the cost function exploits global texture cues in that the contours are associated with salient edges around object boundaries. On the other hand, we look for multiple curvilinear structures, which are not necessarily closed, are latent in the homogeneous texture. Compared with the contour, the curvilinear structures are estimated by subtle local image features. Internal similarity of the structure and an accurate design of shape prior are essential to solve our problem.

Contributions and outline

In this dissertation we have an ambitious goal to propose a general framework which can analyze varied types of curvilinear structures. We assume that the curvilinear structure can be decomposed into a set of line segments. We hence look for an optimal set of line segments which corresponds to the latent curvilinear structure. We mainly proposed two curvilinear structure reconstruction algorithms based on stochastic model with sampling technique (Chapter 3) and ranking learning system with graph representation (Chapter 4). For the stochastic model, we propose a posterior probability of line segments for a given image, and optimize the proposed distribution via Markov chain Monte Carlo (MCMC) sampling. Our work specially contributes to reduce hyperparameter dependency for stochastic process. For this purpose, we simultaneous simulate multiple Markov chains with different parameter settings and aggregate multiple line hypotheses. Moreover, we learn a ranking function which measures the comparability of the given line segments and the latent curvilinear structures. A novel graph representation method using the output rankings of the line segments is proposed to show the topological importance of the underlying curvilinear structure. We test various types of datasets, and provide experiment results in each chapter.

Chapter 2: Curvilinear Structure Modeling Reviews

In Chapter 2, we provide the curvilinear feature descriptors based on image gradient and morphological profiles. 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 comparing to its surroundings, and shows thin and elongated shape. Image gradient is computed with respect to various orientations to detect locally oriented features. Moreover, we employ morphological profiles to find pixels showing dissimilar textures to the backgrounds on an elongated path. In this chapter we also review the previous works on curvilinear structure reconstruction models such as segmentation, graph-based representation, stochastic modeling, and machine learning based approaches.

Chapter 3: Stochastic Model

In Chapter 3, we study a stochastic model for curvilinear structure reconstruction. We propose a probability density of line segments for given image, where this probability density consists of data likelihood and prior term. For data likelihood, we exploit rotated image gradient and morphological filtering responses to estimate the presence of curvilinear structure. We approximate the statistics of curvilinear features using Gaussian distribution with support vector machine. We also constrain local geometric behaviors of line segments. We encourage the smooth connection of line segments, whereas we reject the congestion of line segments. To optimize the proposed density, we employ Reversible jump Markov chain Monte Carlo sampler with delayed rejec-tion. Main contribution of this chapter is to overcome the performance dependency of hyperparameter selection using stochastic optimization. Multiple Markov chains are simulated with different hyperparameter vectors to obtain different line hypotheses, and then we aggregate these hypotheses to reduce sampling space. We provide experiments to compare the curvilinear structure reconstruction results obtained by the proposed algorithm and competing algorithms.

Chapter 4: Inference of Curvilinear Structure

Chapter 4 discusses the structure inference model which quantifies a score according to input-and-output pair information. For classification, the previous works globally set a threshold to remove pixels which show a low probability to be the curvilinear structure. During this process, most of previous works implicitly ignore the correlated information of the pixels on curvilinear structure. We aim to learn a ranking function of the line segments for the given curvilinear structures. The output ranking is able to arrange the line segments in the order of correspondence to the latent curvilinear structures. Specifically, Structured Support Vector Machine (SSVM) is employed to infer the output rankings. We then visualize the topological orders of the curvilinear structure. We study graph theory to first reconstruct curvilinear path which connects the remotest pixels on the latent curvilinear structure. The proposed algorithm progressively adds fine branches into the main structure until it finds the minimum length of structure. We demonstrate experimental results on various datasets to show the efficiency of the proposed algorithm.

Conclusions and Perspectives

We summarize the work and contributions of this thesis. We also provide future directions of the research.

Appendix A: Datasets

Datasets to test the proposed algorithm and competing methods have been included in the appendix.

Appendix B: Thesis overview in French

For French readers, we provide thesis overview of each chapter in French.

List of Publications

Some of the works described in this dissertation has been published in conference proceedings and submitted to an international journal (IJCV, under review). Please refer to the following list of publications. T HE existing curvilinear structure reconstruction models are employed in many different areas, such as:

Submitted Journal

Peer-Reviewed Conference Proceedings

Paper in Preparation

• Medical applications for vascular network analysis;

• Biological image analysis;

• Facial wrinkle detection;

• Satellite image analysis for road network extraction; and Brown 11]. Specifically, graph-based methods have been proposed to represent the curvilinear structure using tree-like structure [Fischler 81, González 10, Türetken 11,Zhao 11, Wang 11,Peng 11,Türetken 13b]. These algorithms first detect a set of seed points corresponding to the local maxima of tubularity measure of the given image. Then, the seed points are connected to reconstruct underlying curvilinear structures. Surplus paths are corrected by pruning wrong edges which induce topological errors [ Peng 11].

Facial wrinkles are an important feature to understand aging progress. In computer vision, facial wrinkles have been studied for age estimation [Fu 10], aged face synthesis [Suo 07, Fu 10], skin texture modeling [ Cula 05] and forensic application [Batool 13]. Facial wrinkle detection algorithms could be involved with the above mentioned applications to analyze aged faces in depth. Localization of facial wrinkles is challenging due to complex patterns of wrinkles and skin textures. In [Batool 12,Jeong 14], authors developed stochastic image models which reconstruct the curvilinear structure (facial wrinkles) with a set of disjoint line segments. They employ the image gradient magnitude to extract linearity feature and constrain the spatial distribution of the line segments.

In remote sensing, road network extraction algorithms have been explored in order to interpret the geographical information [ Defects on road surfaces (road cracks) appear on the thin lines heading in similar direction. For the purpose of convenience and safety, it is preferred to evaluate road surface conditions based on image processing techniques. The task to detect road crack on the images is difficult due to the optical noises and the rough surface textures. Morphological operator to detect cracks is studied with respect to its spatial patterns [ Iyer 05]. Chambon et al. formulated an image segmentation problem following the Markov Random Field (MRF) property to detect irregular shape of the road crack [Chambon 10]. In addition, the acquisition systems which combine optical sensor and laser scanner have been developed to measure the quality of road surfaces. An intensive review of both acquisition systems and image processing techniques for road cracks evaluation can be found in [ Chambon 11],

In the following sections, we first explain the way to compute curvilinear feature descriptors based on image gradient and image morphological profiles (Section 2.1). We then review the existing curvilinear structure reconstruction methods which are segmentation, graph-based model, stochastic model, and machine learning approaches (Section 2.2).

Curvilinear Feature Descriptors

An image I : R 2 → R can be regarded as a function which maps a pixel coordinate into an intensity value. Suppose that an image contains a curvilinear structure of an arbitrary shape. In the previous works, image gradient and morphological features have been employed to describe curvilinear structures which are latent in natural images. The geometric characteristics of a curvilinear structure have important roles and must be taken into account when developing a curvilinear feature descriptor. Geometrically speaking, curvilinear structures are rotatable and have an elongated tube-like shape.

Image gradient information obtained by convolution operation is useful to estimate local orientation and to describe the shape of the curvilinear structure. Also, the curvilinear structures are consecutively distributed toward a similar direction. Since morphological operation is able to recursively detect structural elements of the image, it highlights a sequence of pixels on the underlying curvilinear structure. In the remainder of this section, we review how to implement curvilinear feature descriptors based on image gradient and morphological profiles in detail.

Image gradient profiles

Convolution filtering is a basic operation in image processing to extract low-level image features. For example, Gabor filter banks are popularly employed to detect object boundaries [Mehrotra 92] and to classify fine-grained textures [Jain 90,Grigorescu 02]. Intensive studies to analyze scale space analysis for edge and ridge detection have been done by [ Lindeberg 98]. Specifically, researchers have designed rotatable filters (also known as steerable filters) to estimate local orientation of curvilinear structures [Freeman 91, Perona 95, Jacob 04]. An efficient design of rotatable filters is proposed to reconstruct an arbitrarily oriented kernel as a linear combination of basis kernels. Thus, we perform a set of convolution filters to analyze the given image in different orientations (see Figure 2.1).

The main idea behind such convolution filter design is to create line shaped templates which are analogous to the underlying curvilinear structures. Second order derivative of anisotropic Gaussian kernel is a good approximation to describe shape of the tubular structures. Moreover, the orientation of the anisotropic Gaussian kernel is easily controllable; this property is helpful to estimate the local orientation of the curvilinear structure.

We first analytically develop convolution filter banks to extract curvilinear features.

Let G : [-W 2 , W 2 ] 2 → R be an isotropic Gaussian kernel: G σ 2 : G(x, y; σ 2 ) = exp - x 2 + y 2 σ 2 , (2.1)
where x and y denote the distance from origin (0, 0) in horizontal and vertical axis, respectively. W denotes size of the kernel, and σ denotes the smoothing factor of the Gaussian kernel. Adding an elongation parameter γ > 1 makes the isotropic Gaussian become anisotropic for x axis:

G σ 2 ,γ : G(x, y; σ 2 , γ) = exp - x 2 + (y/γ) 2 σ 2 .
(2.2) The local orientation is taken into account by algebraically performing a rotation in R 2 . That is x θ = x cos(θ) -y sin(θ) and y θ = x sin(θ) + y cos(θ):

G σ 2 ,γ,θ : G(x, y; σ 2 , γ, θ) = exp - x 2 θ + (y θ /γ) 2 σ 2 . (2.3)
Finally, partial derivative of (2.3) with respect to x θ is performed to build a candidate convolution filter

F θ : [-W 2 , W 2 ] 2 →
R that accentuates image gradient magnitudes for θ direction:

F θ = ∂ 2 ∂x 2 θ G σ 2 ,γ,θ (2.4) ∂ 2 ∂x 2 θ G(x, y; σ 2 , γ, θ) = 4x 2 θ σ 4 - 2 σ 2 exp - x 2 θ + (y θ /γ) 2 σ 2 . (2.5)
Alternatively, we show a constructive method to build a rotatable convolution filter banks as the combination of shifted anisotropic Gaussian kernels. An anisotropic Gaussian kernel is shifted to the aperture size a with respect to the minor axis direction as follows:

G σ 2 ,γ,θ,a : G(x, y; σ 2 , γ, θ, a) = exp - (x θ -a) 2 + (y θ /γ) 2 σ 2 . (2.6)
Recall that a candidate filter F θ is Mexican-hat like shape to extract tubular structures. We combine shifted anisotropic Gaussian kernels to imitate the second derivative operation such as:

F θ = G σ 2 ,γ,θ,a + G σ 2 ,γ,θ,-a -G σ 2 ,γ,θ .
(2.7)

Figure 2.2 visually compares a candidate filter created by constructive approach and analytic approach. Comparing these two approaches to create convolution filter banks, the constructive approach (2.7) is apt to estimate the size of curvilinear structures. While the smoothing factor σ of Gaussian kernel is related to the size of the curvilinear structure, it is more likely to blur image signals. On the other hand, the large aperture value is used to highlight pixels on wide curvilinear structure and a small value corresponds to thin structures (see Figure 2.3).

Image morphological profiles

Because image convolution is defined to exploit spatial neighboring pixels, the output is apt to detect local geometric features of the curvilinear structures. However, in the case of uneven illuminations and heavy noises, local image features fail to classify pixels on the same curvilinear structure if they are located far from each other. To make up the disadvantages using local image features to describe curvilinear structure, we compute a morphological operator to exclude pixels corresponding to noncurvilinear structures. To detect thin and elongated shape structure, morphological operations have been designed for medical [Zana 01], remote sensing [Valero 10] and biological [Zhao 11] applications.

The morphological operation is performed to manipulate binary images accord- ing to the structural similarity. Recall that opening, a basic morphological operation, removes objects, which are brighter than their surrounding and smaller than a structuring element. Similarly, path opening operator [ Talbot 07] excludes pixels which do not belong to the curvilinear structures. Let represent an image as a graph, where each pixel corresponds to a node on the graph. The structural elements of curvilinear structure is defined on adjacency graphs with the consideration of topological characteristics, i.e., thin and elongated path. Figure 2.4 (a) shows an adjacency graph consisting of west-east cones (W-E graph) which is used to find horizontal paths. To detect vertical paths, W-E graph is rotated in 90 • and then path opening is applied.

For a given binary image B : R 2 → {0, 1}, path opening seeks for valid paths on the adjacency graph. Let ρ L (u) = {v ∈ B, u → v} be the structuring function which defines a path with a given length L, where u → v denotes a directed edge going from node u to v on the adjacency graph. Then, a valid ρ L -path of length L is represented as L-tuple u = (u 1 , . . . , u L ) that satisfies u k+1 ∈ ρ L (u k ) for k = 1, . . . , L -1. The path opening operator is defined as the union of all valid ρ L -paths of L in the binary image.

To generalize the binary path opening operator to the case of gray-scale images, we decompose an input gray level image into binary images with decreasing order of threshold values. For each decomposed binary image, the binary path opening operator is performed to find the valid paths. The infimum threshold value employed to create the binary image becomes the output of gray level morphological filtering. Finally, morphological profiles are composed of path opening results with different length parameters L. Figure 2.5 shows morphological filtering results on a DNA image.

The curvilinear structure descriptors are employed to measure the quality of reconstructed curvilinear structures based on the image data. However, it is ambiguous to discern pixels corresponding to the curvilinear structure from object boundaries or high frequency noise components. To improve performances, geometric prior and machine learning techniques are involved to advanced curvilinear structure reconstruction modeling. In the following section, we review the state-of-the-art works on curvilinear structure reconstruction modeling.

Existing Curvilinear Structure Reconstruction Methods

To compute a likelihood probability whether a pixel is on the latent curvilinear structure, segmentation techniques can be applied. Low-level image features,e.g., image gradients, are mainly used to estimate this likelihood. On the other hand, graphbased models interpret the curvilinear structure as a sequence of nodes on a graph. They exploit low-level curvilinear features to localize nodes on the curvilinear structure, and define geometric constraints to reconstruct the curvilinear structure with plausible shape. Bottom-up approach represents the curvilinear structure as a disjoint set of line segments. For this purpose, stochastic models have been proposed to define the underlying distribution of the line segments, where an optimal solution consists of smoothly connected line segments. Furthermore, machine learning based approaches are able to detect curvilinear structures on various datasets. We have grouped the existing curvilinear structure reconstruction methods into four classes, which are reviewed in the following sections:

Curvilinear structure segmentation

Curvilinear structure segmentation algorithms can be used to describe irregular shape of the curvilinear structure by evaluating a local tubularity score S(x i ) at pixel x i , i.e., a score corresponds to the likelihood for their pixel to belong to a curvilinear structure. Then, a threshold value is set to decompose the curvilinear structure. To consider varying widths of the curvilinear structure, researchers [Frangi 98, Law 08, Türetken 13a] evaluated their tubularity scores over different scale (width) parameters. The multiscale version tubularity scores S MS (x i ) are maximized at a scale factor r corresponding to the underlying width of the structure at each pixel:

S MS (x i ) = max r min ≤r≤rmax S(x i ; r). (2.8)
We briefly review popular curvilinear structure segmentation algorithms which are mainly applied to vascular network extraction: Hessian matrix analysis [Frangi 98] and Optimally Oriented Flux [ Law 08].

Hessian matrix analysis

In many engineering problems, Taylor approximation is widely employed to analyze the given function as the sum of an infinite series. To interpret the local image structures, Taylor expansion of the image function is given by:

I(x i + δx i ; r) I(x i ; r) + δx i ∇I(x i ; r) + δx i H(x i ; r)δx i ,
(2.9)

where x i = (x i , y i ) denotes the i-th pixel of the given image I, δx i is an infinitesimal distance from x i , r is a scale factor, and ∇I(x i ; r) is the gradient vector of I(x i ; r). Note that image gradient is computed by convolution with derivative of Gaussian kernel. That is

∇I(x i ; r) = ∂ ∂x I(x i ; r) ∂ ∂y I(x i ; r) =   ∂ ∂x G |σ 2 =r * I (x i ) ∂ ∂y G |σ 2 =r * I (x i )   =   ∂ ∂x G |σ 2 =r * I(x i ) ∂ ∂y G |σ 2 =r * I(x i )   .
(2.10)

The Hessian matrix H(x i ; r) is related to the second order derivative of the image, which is computed by

H(x i ; r) =   ∂ 2 ∂x 2 I(x i ; r) ∂ 2 ∂x∂y I(x i ; r) ∂ 2 ∂x∂y I(x i ; r) ∂ 2 ∂y 2 I(x i ; r)   . (2.11)
Intuitively, the second derivative of the image (Hessian) reflects the local image geometries. Hence, researchers examined the eigenvalues of the Hessian matrix to enhance curvilinear structure in the vascular images [Koller 95, Sato 97, Frangi 98]. More precisely, for 2D images, the eigenvalues of the Hessian matrix, λ 1 and λ 2 , show the principle directions of the local structure. With the consideration of the line-shaped structure, pixels on the curvilinear structure are locally aligned toward the same direction. Frangi et al. [Frangi 98] proposed a vesselness measure using the geometric ratio of the eigenvalues R = λ 1 λ 2 , where |λ 1 | ≤ |λ 2 |. Also, the presence of the curvilinear structure (blood vessel) is obtained by the magnitude of the eigenvalues, which are associated with Frobenius norm of the Hessian matrix,

H F = 2 j=1 λ 2 j .
In [Frangi 98], the likelihood of being the curvilinear structure is proposed by the combination of such two observations:

S H (x i ; r) =    0 if λ 2 > 0, exp -R(x i ;r) 2 2β 2 1 -exp - H(x i ;r)| 2 F 2c 2
otherwise, (2.12)

where β and c control the sensitivity of the corresponding measures R and H F , respectively. To take into account different scales (widths) of the curvilinear structure, the vesselness measure (eq. 2.16) introduces the scale factor r which is related to the bandwidth of the filters to compute image derivatives. While the approaches to analyze the eigenvalues of the Hessian for evaluating the local tubularity score of pixels are simple, it is hard to generalize them if the image is heavily corrupted by noise. Also, the results often contain ringing artifacts around curvilinear structure boundaries.

Image gradient flux

Other works, such as [Vasilevskiy 01, Bouix 05, Law 08, Türetken 13a], analyze the local image gradient flux for detecting curvilinear structures. For simplicity, it is assumed that pixels on the curvilinear structures are brighter than backgrounds. Then, at the object boundaries, inward gradient flux is minimized (the outgoing flux is maximized in contrary). Recall that the shape of the curvilinear structure is symmetrical for its cross section and the boundaries are in contrast with the background textures. Therefore, the amounts of the image gradient flux within a local spherical region indi-cate the presence of the curvilinear structure. Meanwhile, the direction of the gradient flux is orthogonal to the local structure.

Law and Chung [Law 08] proposed the Optimally Oriented Flux (OOF) which measures projected image gradient flux along an optimal direction at a local spherical region. In [ Schey 05], the outgoing image gradient flux I f is defined as the surface integral of gradient vectors toward the direction of x over a surface S with radius r:

I f (x i ; r, v) = ∂S ((∇I(x i + rn) • v)v) • ndA = v Q(x i ; r)v, (2.13)
where dA is the infinitesimal area on ∂S, and n = [n 1 n 2 ] is the unit normal vector of the surface ∂S. The matrix Q(x i ; r) encodes the amount of gradient flux projected onto the surface S with the normal vector n, which is given by

Q(x i ; r) = ∂S ∂ ∂x I(x i + rn)n 1 dA ∂S ∂ ∂x I(x i + rn)n 2 dA ∂S ∂ ∂y I(x i + rn)n 1 dA ∂S ∂ ∂y I(x i + rn)n 2 dA . (2.14)
In practice, the matrix Q(x i ; r) is implemented with the set of linear filtering responses, where the linear filters are second derivatives of Gaussian kernel along the corresponding axes. The optimal direction of OOF is computed by minimizing the inward flux at the local spherical structure. It is equivalent to the solution of the following constrained optimization problem: max

v v Q(x i ; r)v subject to v v = 1 . (2.15)
With the Lagrange multiplier method, (2.15) is solved as a generalized eigenvalue problem [ Boykov 04]. Finally, Law and Chung [ Law 08] employed a geometric mean of the eigenvalues of matrix Q(x i ; r) to evaluate tubularity score at x i with scale factor r:

S OOF (x i ; r) =    |λ 1 λ 2 | λ 1 ≤ λ 2 < 0, 0 otherwise, (2.16)
When compared with the methods based on the Hessian matrix analysis, computation of the OOF is less dependent on convolution operations to compute the second derivative of the image.

The curvilinear structure segmentation algorithms finally provide a likelihood probability (tubularity score) of being on the curvilinear structure for each pixel. The tubularity score is employed to enhance the visibility of curvilinear structure on the image. However, it is ambiguous to decide the probability bound based on the observation of the likelihood. For instance, a low threshold value reduces the precision to represent the curvilinear structure. To resolve this issue, geometric priors and global energy optimization are included into curvilinear structure reconstruction modeling.

Graph-based models

Topology of the curvilinear structure often refers to a tree-like structure. Graph based models [Fischler 81, González 10, Türetken 11, Zhao 11, Wang 11, Peng 11,Türetken 13b] reconstruct the curvilinear structure with geodesic paths in a graph G = (V, E), where V = {v i } denotes a set of vertices and E = {e ij = (v i , v j )} is a set of edges. For an efficient computation, the vertex set V (also called seed points) is subsampled from the whole pixels of the given image. Topologically speaking, the seed points are willing to be junctions of the curvilinear structure, and the tubularity scores are locally maximized at these points. In practice, thinning operator [Lee 94] is performed with respect to curvilinear segmentation maps or maximum intensity projection images to detect junction points.

Various methods have been proposed to trace an optimal path that connects seed points on the graph. Sargin2011 et al. developed a Hidden Markov Model (HMM) to estimate a curvilinear path [ Sargin 11]. Peng et al. intially reconstructed an over-complete tree via Dijkstra's algorithm [Dijkstra 59] and discarded redundant nodes [ Peng 11]. Global energy optimization problem is formulated to maximize a posteriori probability of the image data and geometric constraints [González 10, Türetken 11, Türetken 13b]:

t * = argmax t∈G p(t|I, G) (2.17) = argmax t∈G p(I, G|t)p(t), (2.18)
where t is a sequence of vertices which are chosen to reconstruct a curvilinear structure. Thus, the first term p(I, G|t) measures the quality of the reconstructed path based on the image data. p(t) is designed to regularize the geometric properties of the tree structure, e.g., connectivity, orientation consistency, and crossover. The solution of (2.18) can be obtained by either Ant Colony Optimization (ACO) [Türetken 11] or integer programming [Türetken 13b]. Although the geometric priors are employed, an additional relaxation term is needed to modify the topological errors (see Figure 2.6).

Stochastic image models

Curvilinear structures can be seen as a combination of small line segments. Stochastic models define a distribution of line segments with geometric priors to detect multiple line segments corresponding to the curvilinear structure [Lacoste 05, Tu 06, Schlecht 07, Batool 12, Predoehl 13, Jeong 14]. Each line segment s encodes its shape information, i.e., length and orientation. Thus, a line segment s is defined as a tuple of the center coordinate x, length , and orientation θ. To optimize the probability density, Reversible Jump Markov chain Monte Carlo sampler [Green 95] has been typically employed. A state of the Markov chain is related to a line configuration. The Markov chain is locally perturbed by sub-transition kernels, and converges at a stationary state which maximizes the probability density (see Chapter 3 for details). However, the stochastic models require heavy formalization to interpret spatial distribution of the line segments. Large number of parameters should be specified to describe the shape of the line segments and to control the relative importance of data and prior energy terms. Stochastic modeling is considered less practical for applying to a large amount of varied datasets because its performance is very sensitive to the selection of parameters.

That is s = (x, , θ) ∈ I × R × Θ,

Machine learning based approches

Previous works on curvilinear structure modelling have been developed independently to tackle a specific application. Recent trends on computer vision research proceed to handle many data as possible. Thus, it is desired that the solution should be self-contained and use the minimum number of user-defined parameters. Apparently, for this purpose, machine learning techniques are involved with many computer vision problems nowadays. Specifically, model parameters can be learned automatically with labeled dataset (training set), so that we can reduce efforts to adapt the proposed model for different types of dataset. Becker et al. [Becker 13] opened the possibility to employ a machine learning technique for curvilinear structure segmentation. They proposed a Gradient boosting framework to obtain optimal convolution filter banks to extract the curvilinear features. Scale space distance transform (SSDT) is proposed as a regression function to predict the scale (width) of the curvilinear structure and localize centerlines [Sironi 14] .

In the general setup of the machine learning problems, we aim to learn a function h : X → Z between an input space x ∈ X and output space z ∈ Z. To evaluate the quality of a prediction ẑ = h(x), a loss function ∆(z, h(x)) is defined. Then, a machine learning algorithm minimizes the loss function to find an underlying distribution of data with a given training dataset

D = {(x i , z i )} N i=1
, where N is the number of input and output pairs.

Gradient boosting [Zheng 07, Hastie 09] approximates a function as a linear combination of weak classifiers:

h(x) = M j=1 α j ϕ j (x),
(2.20)

where x ∈ R n denotes an input vector, α j ∈ R is a weight of the weak classifier ϕ j (•) : R n → R, and M is the number of week classifiers. For a given training dataset D, an optimal parameter set of regression tree h is obtained in a greedy manner.

In other words, at each weak classifier ϕ j , the corresponding weights are found via exhaustive search on a regression tree to minimize a total loss L = N i=1 ∆(z i , h(x i )).

Supervised curvilinear feature learning

In [Becker 13], Gradient boosting framework is employed to create convolution filter candidates with a closed form, and to learn a regression tree for the curvilinear structure segmentation. For the curvlinear feature learning and classification, Gradient boosting implementation of [Becker 13] consists of two phases: (1) searching for a set of convolution kernels and (2) learning parameters of the regression tree which minimize the splitting cost. A weak classifier ϕ j (x i ) corresponds to the inner product of a convolution kernel Fp and a vectorized image patch X i at a sampled point x i . Thus, a candidate kernel is found:

Fp = argmin F N i=1 ω i (F X i -g i ) 2 + m∼n F (m) -F (n) 2 , (2.21)
where ω i and g i denote weight and target values, respectively. The partial differentiation of the loss function is computed to obtain ω i and g i . The second term of eq. (2.21) induces smooth kernel values with local neighborhoods m ∼ n, where F (m) denotes m-the value in the filter F. Then, regression parameters for decision tree are computed with respect to each Fp .

Centerline detection based on learning a scale-space distance transfrom

Gradient Boosting framework is exploited in a different way to localize a centerline of the curvilinear structure in [Sironi 14]. In this work, the centerline also encodes the scale of the curvilinear structures. Previous curvilinear structure segmentation methods often highlight pixels more on the boundaries of the structures than on the centerlines, even though the models take into account multiscale terms. Thus, in [Sironi 14], the authors proposed a regression model which learns the scale of the curvilinear structure based on the distance transform.

Let D C (x i ) denote the distance transform which measures the distance between the given point x i and the closest point belonging to the centerline. To train a regressor, the distance transform D C (x i ) is updated into a function d(x i ), whose values monotonically decrease from the centerline:

d(x i ) =    exp 6 1 -D C (x i ) d M -1, if D C (x i ) < d M 0 otherwise, (2.22)
where d M truncates the values according to the size of convolution filters. The configurations for learning a regressor are given by

D = {(f i , z i )} N i=1
, where f i ∈ R n denotes feature vector of x i on the image I and z i = d(x i ). Specifically, the feature vector concatenates convolution filtering responses

f i = [(F 1 * I) (x i ), • • • , (F n * I) (x i )] ,
where the convolution filter banks are computed by [Rigamonti 12]. Gradient Boosting [Hastie 09] is used to learn a regression function. Finally, likewise Canny edge detector [ Canny 86], Non Maximum Suppression (NMS) is performed along to the perpendicular direction of the structure for thinning the output function.

Conclusions

The previous works on curvilinear structure reconstruction modeling have been studied geometric features of the curvilinear structure on images. We review the curvilinear structure reconstruction algorithms to analyze blood vessel networks, road net-works, biological images, facial wrinkles, and defects of materials. To summarize, the curvilinear structure locally exhibits thin, elongated, and tube-like structure on images. We summarize the implementation of curvilinear feature descriptors based on image gradient profiles and morphological profiles. Image gradient information is used to extract local image features which are analogous to elongated tubular structures and to estimate local orientation. A set of convolution filters is proposed to capture such geometric properties. Morphological profiles detect valid paths on a graph which reflects the topological feature of the curvilinear structure. We performs the morphological operation to find remote pixels on the same curvilinear structure within a pre-defined length.

Topology of the curvilinear structure varies depending on its application, so that different representation techniques are preferred to reconstruct the underlying curvilinear structures. Segmentation algorithms are apt to analyze irregular shaped curvilinear structure and provides a tubularity score for each pixel. Graph-based algorithms represent the latent curvilinear structure as a geodesic path on a graph. Stochastic image modeling reconstructs curvilinear structures as a set of disjoint line segments. Geometric prior defines the spatial interaction of the line segments. Machine learning based approaches is developed to analyze curvilinear structure not depending on a specific application.

In this dissertation, we aim to develop curvilinear structure reconstruction methods which are applicable to various types of datasets. In Chapter 3, we propose a stochastic model which employs image gradient profiles as the image based term and constrains spatial interaction of line segments. In Chapter 4, we develop a ranking learning system and graph based representation to infer curvilinear structures.

Chapter 3 Stochastic Modeling of Curvilinear Structure Reconstruction Abstract

This chapter presents a novel stochastic model and the corresponding optimization technique to reconstruct curvilinear structures. We sample multiple line segments which are associated with the latent curvilinear structure. Therefore, we aim to maximize a posterior probability density with respect to a set of line segments for given image data. To fit the set of line segments into image data, we compute curvilinear feature descriptors based on oriented image gradient and morphological profiles. We employ support vector machine to estimate the distribution of the proposed curvilinear features and to regularize sampling space. We constrain local interaction of the line segments to obtain smoothly connected line configuration. The optimization technique consists of two steps to reduce the significance of the parameter selection in our stochastic model. We simulate several Markov chains with different parameter settings to collect line hypotheses on the same configuration space. Then, we maximize the consensus among line hypotheses to reduce the sampling space and to improve the reliability of the curvilinear structure reconstruction. Experimental results show that the algorithm effectively localizes curvilinear structures on a wide range of datasets. Stoyan 87, Lieshout 00, Descombes 02, Chatelain 12] is a stochastic process which can enforce high level constraints on shape prior. However, the MPP model requires heavy formalization to interpret spatial distribution of the objects. Large number of parameters should be defined to describe the geometric shape of the objects (modeling parameters) and to control the relative importance of data and prior energy terms (hyperparameters). MPP modeling has been considered less practical to solve general problem because the performance is very sensitive for the selection of parameters. Although stochastic expectation maximization algorithm [Celeux 96, Møller 03] has been used to estimate modeling parameters, it exhibits both speed and scalability issues.

In this chapter, we propose a novel stochastic model based on MPP framework for curvilinear structure extraction in a fully automatic way, where the performance is not biased by the hyperparameter selection. Indeed, our stochastic model can detect wide types of input data without a sophisticated parameter tuning process. In Section 3.1, we review MPP framework [Stoyan 87, Lieshout 00] to provide a mathematical description of the proposed model. To fit into image data, we proposed curvilinear features with rotated image gradient and morphological profiles. Support vector machine (SVM) [Vapnik 95] is employed to approximate the data likelihood distribution of the given proposed curvilinear features. The prior distribution is also defined on local configuration to implement smooth connection among line segments (Section 3.2). Markov chain Monte Carlo (MCMC) sampler [Gilks 95, Green 95, Robert 04] with delayed rejection scheme [Green 01] is employed to optimize the proposed probability density function. We explain the implementation details of transition kernels for MCMC sampler in Section 3.3. To avoid the burden of hyperparameter selection, we generate multiple candidates of the line configuration with different hyperparameter settings. Next, we combine the whole set of line candidates in a way that maximizes the consensus among line detection results (Section 3.4). Extensive experiments on various datasets including facial wrinkles, road cracks, DNA filaments, and blood vessels demonstrate the effectiveness of the proposed model for extracting thin curvilinear structures (Section 3.5).

Revisited Marked Point Process Framework

Modeling with spatial point process aims to describe the spatial pattern of data using a finite number of parameters. We evaluate the spatial point patterns with respect to randomness in a finite space. For example, we employ a statistical model to realize a completely independent point pattern. Then, we observe how the given points are differently (or similarly) distributed in space comparing to the randomly generated point patterns (see Figure 3.1).

Formally speaking, a spatial point process X is a random variable and its realization is associated with an unordered set of points Ω n = {x 1 , . . . , x n } in a compact set F ⊂ R d in d-dimensional space. A spatial point process on F maps from a measurable space (F, B, ν) onto the set of configurations Ω = ∪ ∞ n=0 Ω n , where B denotes σ-algebra of B, and ν is a positive counting measure on F. In other words, for all bounded Borel set B ⊆ F, the number of points N X (B) falling in B is a finite random variable.

Furthermore, for each point we add an additional information to describe the shape of an object. Let s i = (x i , m i ) be an augmented point (object) which encodes its center location x i ∈ F and a shape information (mark) m i ∈ M. Similarly to sample a point with probability measure µ, we obtain a parameter to describe its shape with a probability measure µ M . Thus, a realization of the marked point process on Ψ = F ×M is a finite random configuration set s = {s 1 , . . . , s n } ∈ Ψ.

In practice, the homogeneous Poisson process is employed to develop a stochastic model to quantify randomly scattered point patterns. The homogeneous Poisson process X with a positive intensity ν > 0 satisfies the following properties: Property 1. For all bounded Borel sets B, the number of points N X (B) falling in B follows a discrete Poisson distribution with mean ν|B|:

Pr(N X (B) = k) = ν|B| k k! e -ν|B| .
Property 2. For disjoint bounded sets, B 1 , . . . , B k , the corresponding random variables N X (B 1 ), . . . , N X (B k ) are independent.

Straightforwardly, the Poisson process can be generalized by a non-constant intensity measure ν(•) on F. This process is called an inhomogeneous Poisson process. The Poisson processes realize independently scattered points on the space, while real objects can be interactive on a local configuration. For instance, we want to avoid overlapping between objects and to regularize the orientations. To take into account interactions of the scattered points, we specify a density f (•) with respect to the homogeneous Poisson process with intensity measure ν(•). Let f (•) be a non-negative function on sampling space. Then, the measure ν(•) having density f (•) with respect to ν(•) is defined by:

ν(B) = B f (s)ν(ds). (3.1)
For example, Markov point process defines f (•) with the pairwise interaction (repulsion or attraction) of points according to their distance. Figure 3.2 compares the realization of homogeneous Poisson process, inhomogeneous Point process, and repulsive Markov point process, respectively. Now we design the density function f (•) which can control point pattern according to the given data and prior knowledge of their physical properties. The density function also refers to a probability distribution of the points. In this chapter, we represent a line segment s i = (x i , i , θ i ) as a tuple of its center location, length, and orientation, respectively. To reconstruct the curvilinear structure, we collect smoothly connected line segments on the configuration space Ψ. In other words, our goal is to find an optimal set of line segments ŝ which maximizes posterior probability for given image f (s) p(s|I). The posterior probability can be separated into data likelihood p(I|s) and prior probability p(s) such as:

ŝ = argmax s∈Ψ p(s|I) = argmax s∈Ψ p(I|s)p(s). (3.2)
In general, Markov chain Monte Carlo (MCMC) samplers [Green 95, Green 01, Gamal-Eldin 10, Verdié 14] are employed to maximize the proposed density function p(s|I). Each state of a discrete Markov chain (X t ) t∈N corresponds to a random configuration on the Ψ. The chain is locally perturbed by transition kernels, and is evolved to converge to the stationary distribution which is identical to the maximum of the proposed probability density. In Section 3.3, we further explain MCMC dynamics and implementation of the transition kernels.

Stochastic modeling of curvilinear structure

In this section, we propose a stochastic model for curvilinear structure reconstruction. The stochastic model consists of data likelihood p(I|s) to localize a set of the line segments s for given image data I and prior distribution p(s) to define local geometry of line segments. Let z i ∈ {-1, 1} denote the latent variable for pixel x i ; z i = 1 if x i is on the curvilinear structure, otherwise z i = -1. We predict a pixel x i whether it belongs to the curvilinear structure using the corresponding curvilinear feature vector φ(x i ). Specifically, we propose the curvilinear feature descriptors φ : I → R N based on image gradient or morphological profiles, where N represents the dimension of the feature space. The data likelihood combines information of pixels on the given line segment:

p(I|s i ) = 1 Z exp -E data (s i ) = 1 Z exp   - x j ∈s i E data (φ(x j ))   , (3.3)
where E data denote Gibbs energy form and Z = s∈s exp -E data (s) is a normalizing constant. To compute prior probability, we assume that an interaction of the line segments is valid if and only if they are close enough in the image coordinate. Thus, we consider pairs of line segments such that their center distance is smaller than half the sum of their lengths:

s m ∼ s n = (s m , s n ) ∈ Ψ 2 : 0 < x m -x n 2 ≤ m + n 2 + , ( 3.4)
where denotes the marginal distance to be connected with each other. With the consideration of neighboring condition (3.4), we can compute the prior probability as follows:

p(s) = p(s 1 , . . . , s #(s) ) = sm∼sn p(s m |s n ), (3.5)
where #(s) denotes the number of line segments in the configuration s.

The probability density is interchangeable with Gibbs energy form as a negative log probability [Gibbs 02]: E(x) ∝ -log p(x). We denote energy form of data likelihood probability as E data and prior probability as E prior . Therefore, maximizing the proposed posterior probability is equivalent to minimizing the total energy of line configuration:

ŝ = argmax s∈Ψ p(s|I) = argmax s∈Ψ #(s) i=1 p(s i |I) (3.6) = argmax s∈Ψ #(s) i=1 p(I|s i ) sm∼sn p(s m |s n ) (3.7) = argmin s∈Ψ #(s) i=1 E data (s i ) + sm∼sn E prior (s m , s n ), (3.8)

Data likelihood

Data likelihood employs a curvilinear feature descriptor to measure the fidelity of the pixel with respect to the given image data. First, we exploit rotated image gradient magnitudes. A convolution filter F θ is applied to extract tubular structures oriented to the corresponding direction θ. The rotated image gradient information is suitable to detect local tubular structure of the image; however, image gradient often fails to classify pixels even though they are remotely connected in the same curvilinear structure. To make up for such failure case, we propose a curvilinear feature descriptor using morphological profiles. We then combine several feature maps which highlights the latent curvilinear structure with different statistics. Hence, Support vector machine (SVM) [Vapnik 95] is employed to obtain relative importance of each feature map and to approximate the data likelihood distribution of the given curvilinear feature vector.

Curvilinear feature based on image gradient statistics

We perceive the thin linear shape structure due to a contrast dissimilarity with the background textures. Specifically, pixels on the curvilinear structure are slightly darker or brighter than its surroundings, and spread on a similar direction. Thus, we analyze image gradient since it can measure local intensity variation. We compute the rotated gradient information by convolving the given input image I with rotatable filters {F θ }, which are defined in the previous chapter (Section 2.1). Let ∇ θ I = F θ * I denote the image gradient map which accentuates image gradient magnitude corresponding to the orientation θ. When the input image is heavily corrupted by noise or composed of uneven textures, observing gradient distribution often fails to detect linear structures.

To ease this problem, we also assume that the intensity on the line segments is homogeneous. We then compute a curvilinear feature vector φ 1 (x) for pixel x by combining these information :

φ 1 (x) = [I(x), ∇ θ I(x)] ,
(3.9)

Curvilinear feature based on morphological profiles

Inner similarity of pixels on the curvilinear structure is hard to get using local gradient features due to irregular illuminate condition, rough background texture, and noise. Morphological operator detects linear shape features using recursive structural patterns on the image. Hence, we employ path opening operations [Talbot 07] by controlling the length parameter to obtain different feature maps which can compensate the drawbacks of image gradient features. We assume that pixels belonging to the curvilinear structure are brighter than their neighborhoods. Otherwise, we invert the intensity level of the input image before applying the morphological operator. We then denote the curvilinear feature vector φ 2 (x) which combines the morphological filtering responses I L with increasing order of length parameter L. That is

φ 2 (x) = I L 1 (x), . . . , I L 4 (x) , where L 1 < • • • < L 4 .

Initial classification of the latent curvilinear structure

We concatenate two curvilinear feature vectors φ 1 (x) and φ 2 (x) to compute data likelihood of being curvilinear structure:

φ(x) = φ 1 (x) φ 2 (x)
. In this chapter we employ the SVM for the initial classification of the latent curvilinear structure. Let D = {(φ(x i ), z i )} K i=1 be a training set. In practice, we randomly select 50 training data points for positive and negative samples, respectively. To obtain a linear hyperplane, we solve following quadratic programming problem: The solution of (3.10) w defines the decision function h i = h(φ(x i )) = w φ(x i ) + b of the SVM classifier, and estimate the latent label ẑi = sign(h i ). The distance from the decision boundary (hyperplane) to a feature vector can be considered as a psuedoprobability belonging to curvilinear structure, which can be approximated by a sigmoid function [Platt 00]:

min b,w,ξ 1 2 w w + C K K i=1 ξ i subject to z i (w φ(x i ) + b) ≥ 1 -ξ i ; for i = 1, . . . , K (3.
p w (z i = 1|φ(x i )) = 1 1 + exp(αh i + β) , ( 3.11)
where regularization parameters (α and β) are estimated by maximum likelihood. We want to regulate that data likelihood probability is likely to the initial estimation obtained via SVM classifier. Moreover, we assume that the data likelihood of the given curvilinear feature vectors follows Gaussian distribution with mean µ d and covariance Σ d :

µ d = i p(z i = 1|φ(x i ))φ(x i ) i p(z i = 1|φ(x i )) , Σ d = ww 1 2 .
(3.12)

We then define the data likelihood for given curvilinear feature vector φ(x):

E data (φ(x i )) = -ω d log p w (z i = 1|φ(x i )) -log N (φ(x i ); µ d , Σ d ), (3.13) ω d exp(αh i + β) + φ(x i ) -µ d Σ d -1 φ(x i ) -µ d (3.14)
where ω d denotes a hyperparameter to support the initial classification results of SVM. The binary SVM classifier fails to classify the curvilinear structure so that it yields jittering artifacts around boundaries. This is because the SVM classifier assumes identical-independent-distributed (i.i.d.) condition for the sample points when it learns hypothesis. To improve the classification performances, we employ prior constraints in a local line configuration (see Figure 3.4).

Prior Energy

In this section, we propose the prior energy to define spatial interactions of a local line configuration. We define the coupling feature vector c mn for adjacent line segments s m ∼ s n by neighboring condition (3.4). The coupling feature vector takes into account the connectivity and curvature of pairwise line segments, such as:

c mn = [d mn , θ mn ] ,
(3.15)

where d mn denotes the minimum distance between endpoints of s m and s n , which is related to the connectivity. θ mn is the angle difference to measure the curvature of two line segments (see Figure 3.5). Intuitively, small values of d mn and θ mn are desired to obtain connected and aligned line configuration: c mn c mn 0. We also compute the overlapping area Υ(s m , s n ) to reject congestion of lines within a local area. In order to evaluate an overlapping area between line segments, we dilate the line segments with a three pixel-radius disk, and then count the number of pixels falling in the same area. Suppose that we have a set of points A(s m ) which is a dilated version of the line segment s m , and |A(s m )| denotes the number of pixels in A(s m ). As shown in Figure 3.6, we penalize a configuration {s m , s n }, when a portion of the overlapping area is greater than 10% of min{|A(s i )|, |A(s j )|}. The criteria for rejection are then given as Assuming that c mn follows a mixture of Gaussian distributions with zero mean 0 and covariance Σ p , we can write the prior probability:

Υ (s m , s n ) =    0 if |A(s i )∩A(s j )| min{|A(s i )|,|(A(s j )|} < 0.1, 1 otherwise. (3.16) (a) Υ(s m , s n ) = 0 (b) Υ(s m , s n ) = 1
E prior (s m , s n ) = ω p Υ(s m , s n ) - D d=1 π p d log N (c mn ; 0, Σ p d ),
(3.17)

ω p Υ(s m , s n ) + D d=1 π p d c mn Σ p d -1 c mn , (3.18) ω p Υ(s m , s n ) + π 1 d 2 mn + π 2 θ 2 mn , (3.19)
where D = 2 denotes the number of Gaussian mixtures, and π p d is the corresponding weighting coefficient. ω p is a hyperparameter which controls the rejection rate.

The hyperparameters of data likelihood ω d and prior probability {π P d } D d=1 and ω p play important roles to sample reliable line segments. For example, ω d penalizes the birth of a line segment in the area which is not associated with the initial estimation. {π P d } D d=1 promotes adjacent line segments to be connected and aligned in the final configuration. ω p supports the sensitivity of overlapping criteria. However, it is hard to estimate such hyperparameters because of implicit dependencies among them. Specifically, for the prior distribution, we cannot exploit true statistics of related parameters for learning phase. We will tackle this issues in section 3.4.

Optimization

To maximize the proposed probability density of line segments, we employ a Markov chain Monte Carlo (MCMC) sampler. In this section we discuss how the MCMC sampling devotes to maximize the probability density function in a complex sample space and its implementation details.

Figure 3.7: The latent curvilinear structure, denoted by gray color, is reconstructed by a set of line segments. A state of the Markov chain is related to a line configuration, and updates its state according to the transition kernels. We finally obtain an equilibrium state of the Markov chain which maximizes the corresponding probability density in the sampling space.

Markov chain Monte Carlo sampler

We employ a discrete Markov chain {X t } t∈N to support the proposed probability density of line segments in the configuration space Ψ. A state of the Markov chain X t at time t is a random variable which is associate with a set of line segments s t (see Figure 3

p(s a |I)K(s a , s b ) = p(s b |I)K(s b , s a ) (3.21)
The transition kernel consists of several sub-transition kernels which defines different types of jumps to perturb the state:

K(s a , s b ) = m q m K m (s a , s b ), m q m = 1, (3.22)
where q m denotes a probability to choose m-th type of sub-transition kernel K m (s a , s b ).

α m (s a , s b ) = min 1, Q m (s b |s a )p(s b |I) Q m (s a |s b )p(s a |I) . (3.23)
In the following section, we design the proposal kernels which efficiently perturb the line configuration in order to maximize the proposed probability density.

Reversible Jump Markov chain Monte Carlo sampler

Since the dimensionality of the line configuration can be changed during simulation process, we therefore employ Reversible jump Markov chain Monte Carlo (RJMCMC) sampler [Green 95] to maximize the probability density function. The RJMCMC sampler is an iterative method that locally perturbs a current configuration s with a transi-tion kernel, i.e., birth-and-death (BD) and linear transform (LT). Each sub-transition kernel consists of pair of reversible jumps and holds the detailed balance condition (3.21). Acceptance ratio α m (s, s ) is compared with a stochastic value rand[0, 1] to update a new configuration s . The RJMCMC sampler is coupled with the simulated annealing (SA) algorithm [Kirkpatrick 83] to secure the convergence of the Markov chain via relaxation parameter T (temperature); the temperature gradually decreases as the iteration goes on. To compute an acceptance ratio of the transition kernel, we use a density p(s|I) 1 /T instead of p(s|I). The acceptance ratio is

α m (s, s ) = min 1, Q m (s|s ) Q m (s |s) p(s |I) 1 /T p(s|I) 1 /T . (3.24)
Algorithm 1 summarizes the optimization process using RJMCMC.

Algorithm 1 RJMCMC sampler [Green 95]

1: Initialize:

X 0 ← s 0 (or X 0 ← ∅), t ← 0, T ← T 0 2: while T > T min do 3: s ← X t 4:
Choose a transition kernel K m according to probability q m 5:

Propose a new configuration s with Q m (s |s)

6:
if α m (s, s ) > rand[0, 1] then 7:

X t+1 ← s 8: end if 9: t ← t + 1 10:
Decrease the temperature: T ← T t 11: end while More precisely, the BD kernel changes the dimensionality of the current configuration s by adding a new line segment or removing an existing line segment. When the birth kernel proposes a new configuration s = s ∪ {s}, the length and the orientation of the new line segment are uniformly sampled from the mark space M = [ min , max ] × [θ min , θ max ], where min and max are the minimum and maximum length of the line segment, respectively. Similarly, θ min and θ max denote the minimum and maximum orientation of the line segment, respectively. Note that we refuse a birth of the line segment lying on singular points, which have zero gradient magnitudes. On the other hand, the death kernel removes a line segment which is randomly picked from the current configuration. Thus, a new configuration s = s\{s} is proposed by the death kernel. We compute the acceptance ratio of the birth kernel α B and the death kernel α D , given by The LT kernel chooses a line segment s randomly, and then modifies its model parameters: s = (x, ( , θ)) → s = (x ± ∆x, ( ± ∆ , θ ± ∆θ)), where ∆x, ∆ , and ∆θ denote changes of center position, length, and orientation, respectively. The LT kernel draws a new configuration s = s\{s} ∪ {s }. The acceptance ratio of the LT kernel is defined by Figure 3.8 shows how the sub-transition kernels update the corresponding configurations.

α B (s, s ) = min   1, p D p B µ M (M ) #(s) + 1 p (s |I) 1 /T p(s|I) 1 /T   , (3.25) α D (s, s ) = min 1, p B p D #(s) µ M (M) p(s |I) 1 /T p (s|I)
α LT (s, s ) = min 1, p(s |I) 1 /T p(s|I) 1 /T . (3.27)

Delayed Rejection

The LT kernel can be extended by the delayed rejection scheme [Green 01]. In order to reduce the burn-in time, we add heuristics to design the delayed rejection kernel. The main idea of the delayed rejection scheme is to give a second chance to a rejected configuration by enforcing connectivity of the line segments. The acceptance ratio of delayed rejection is defined by

α 2 LT (s, s , s) = min 1, Q LT (s |s) Q LT (s |s) Q 2 LT (s|s , s) Q 2 LT (s|s, s ) [1 -α LT (s, s )] [1 -α LT (s, s )] p(s|I) 1 /T p(s|I) 1 /T , min 1, p(s|I) 1 /T -p(s |I) 1 /T p(s|I) 1 /T -p(s |I) 1 /T . (3.28)
where s = s \{s} ∪ {s }, s = s\{s} ∪ {s}, and Q 2 LT (s|s, s ) is the transition kernel for the delayed rejection. Intuitively, an alternative configuration s is accepted if the corresponding probability p(s|I) is higher than p(s |I). At the same time, it satisfies p(s |I) < p(s|I), which is a good reason to reject the intermediate configuration s . To propose an alternative line segment s, we look for the closest endpoints from both ends of s , which is rejected from the first trial. The line segment s is generated by interpolation of the retrieved points; we force the connectivity of the neighboring segments, so that a probability of being accepted increases in terms of prior distribution. Figure 3.9 summarizes the process of the delayed rejection kernel, and Algorithm 2 provides the pseudo-code of the RJMCMC sampler with delayed rejection.

Algorithm 2 RJMCMC sampler with delayed rejection

1: Initialize:

X 0 ← s 0 (or X 0 ← ∅), t ← 0, T ← T 0 2: while T > T min do 3: s ← X t 4:
Choose a transition kernel K m according to probability q m 5:

Propose a new configuration s with Q m (s |s) Decrease the temperature: T ← T t 18: end while

6: if α m (s, s ) > rand[0,

Integration of Line Hypotheses

While MPP modeling allows to design complex prior knowledge of the object distribution, its performance is very sensitive to the selection of modeling parameters and hyperparameters. For clarity, we note that the modeling parameters are related to the physical characteristics of the line segments (e.g., range of length and orientation). The hyperparameters denote the weighting coefficients of energy terms (i.e., ω d , π p , and ω p ). The modeling parameters can be chosen empirically since the values are related to the image resolution (see Section 3.5); however, it is hard to estimate the hyperparameters via trial-and-error for different types of dataset. Our goal is to maximize the proposed probability density without estimating hyperparameters.

Generation of K Line Hypotheses

Let ω = [ω d ω p π 1 π 2 ] be a hyperparameter vector which consists of the weighting coefficients of the proposed probability density. Suppose that we have

K different hyperparameter vectors, ω 1 , • • • , ω k , • • • , ω K .
For each hyperparameter vector, we substitute k-th hyperparameter vector ω k into the proposed probability density p(s|I, ω k ). Then, we look for its optimal configuration ŝk = argmax s⊂Ψ p(s|I, ω k ) via Monte Carlo sampler proposed in Section 3.3.3.

For the implementation, we bound the values of ω due to the practical reason. More precisely, we fix ω p as 1000 and control the other values according to the proportion of χ = ω d = 1. We set ω 1 = [χ, 1000, -0.1χ, -0.1χ] , and then linearly increase χ to create hyperparameter vectors, i.e., ω 2 = [1.1χ, 1000, -0.11χ, -0.11χ] , ω 3 = [1.2χ, 1000, -0.12χ, -0.12χ] , and so on. In our experiments, we use K = 15 to create line hypotheses.

Algorithm 3 Intergraion of K line hypotheses

1: for all ω 1 , . . . , ω K do 2:

Look for an optimal set of line segments with k-th hyperparameter vector

Combination of Line Hypotheses into a Probability Map

We now have a family of line hypotheses S = K k=1 ŝk obtained from K different hyperparameter vectors. We jointly use the image data and these line hypotheses. More specifically, the final solution s * maximizes not only the probability density w.r.t. image but also the consensus among line hypotheses. For each optimal configuration ŝk , we compute a probability map P k of being a line in the image site. Then, we integrate K probability maps into a mixture density P S :

P k (x) =        1 if ∃s k i ∈ ŝk , x ∈ s k i , 1 2 if ∃s k i ∈ ŝk , x ∈ A(s k i ), 0 otherwise, (3.29) P S (x) = 1 K K k=1 P k (x).
(3.30)

Figure 3.10 compares image gradient magnitude, morphological filtering [Talbot 07], supervised feature learning [Becker 13], line hypotheses, and the mixture density. Since the input image contains many high frequency components, its gradient also highlights non-linear structures in the background. While the morphological filter accentuates linear structures, its performance depends on the setting of path length. Supervised learning method requires high quality of a training dataset and corresponding ground truth images. Depending on the setting of hyperparameter vectors, the MPP model leads incomplete detection results as shown in Figure 3.10. (e)-(g). We integrate line hypotheses of the proposed MPP model into a mixture density P S . The mixture density shows the consensus between line hypotheses in the sense that the pixels corresponding to line structures are more highlighted when compared to [Becker 13, Talbot 07].

We assume that the most promising hyperparameter vector draws a configuration which is more akin to the mixture density. We compute the correlation-coefficient (CC) between P S and P k 's to analyze coherence of line detection results. That is k = argmax k={1,...,K} CC(P S , P k ), (3.31)

CC(P S , P k ) = x (P S (x) -E [P S ]) (P k (x) -E [P k ]) x (P S (x) -E [P S ]) 2 x (P k (x) -E [P k ]) 2 , ( 3.32)
where k represents the index of the most reliable hyperparameter vector.

Curvilinear Structure Extraction from Reduced Sampling Space

The line hypotheses span a configuration space S ⊂ Ψ which will be considered as a new sample space. Since the size of S is significantly reduced compared to the original sample space Ψ, the optimization process becomes more tractable in terms of convergence time and detection accuracy. We redefine the data likelihood energy by adding a new energy term as follows:

E data (φ(x)) = E data (φ(x)) -log P S (x), (3.33) E data (s i ) = x j ∈s i E data (φ(x j )). (3.34)
We simulate the modified probability density over the reduced sample space S with the most promising hyperparameter vector ω k:

s * = argmax s⊂S p(s|I, ω k) (3.35) = argmin s⊂S #(s) i=1 E data (s i ) + sm∼sn E prior (s m , s n ; ω k).
(3.36) We illustrate the line hypotheses integration process in Algorithm 3.

Experiments

We test the proposed algorithm on a wide range of datasets: facial wrinkles, DNA filaments1 , road cracks, and retinas. The facial wrinkle images are collected on the Internet, and forehead areas are manually selected for the experiments. Test images of the defects on the road pavements and ground-truth are courtesy of Chambon et al. [Chambon 10]. We use the DRIVE dataset [ Staal 04] to test the proposed algorithm on retina images.

For all test sequences, we fix the modeling parameters as follows: min is set to 5 pixels and max = 20 pixels. The orientation θ is varying from -90 • to 90 • with increments of 2 • . The marginal distance of connected segments is fixed to 2 pixels, and the maximum angular difference of aligned segments τ is 30 • . For the SA, the initial temperature T 0 is set to 100, and it follows the logarithm cooling schedule T t = T 0 /log(1+t), where t denotes the number of the current iteration. We start the sampling process with the empty configuration. However, careful choice of initial segments can speed up the convergence of the algorithm (see Figure 3.11). The computational time depends on the image resolution and the number of objects; it takes less than a minute for the experimental images having 300 × 400 pixels, approximately. We use a PC with a 2.9 GHz CPU (4 cores) and 8 GB RAM.

To compare the performances of the proposed method with the state-of-the-art techniques, we apply the path opening operator [Talbot 07] on the gradient magnitude images by controlling the length parameters. For the supervised feature learning algorithm [Becker 13], we train 15 images for each dataset. The baseline MPP is selected from the line hypotheses among which it shows the best performance. In our experiments, we use the original implementations of path opening operator2 and supervised feature learning algorithm3 . Figure 3.12 shows the precision-and-recall curves for four test images. To plot the precision-and-recall curve, we need to control threshold values for changing the operating point of the binary classifiers. However, the proposed stochastic algorithm generates a binary map which denotes the areas corresponding to the final line configuration. The performances of the supervised feature learning algorithm [Becker 13] are controlled by the quality of the training set; hence, it shows low performances on wrinkle and dna datasets, which are composed of noisy images of various sizes. In particular, the ground truth set of the wrinkle dataset is based on subjective perception. While the morphology operator [ Talbot 07] enhances linear structures on gradient magnitude images, it is required to specify the length of the linear structures according to the target applications. It might be unfair comparison with the competing algorithm using the precision-and-recall curve. One reason we have to consider is that pixelwise comparison cannot evaluate the structural similarity between ground truth map and the detected line segments. In other words, the measure can mislead the performance if the line detection results are slightly sifted compared to the groundtruth map. For example, although the proposed algorithm provides visually plausible shape reconstruction for road cracks, the quantitative performance is not the best. We visually compare the curvilinear structure detection results of the proposed MPP model with the manually labeled image by human expert, morphology operator [Talbot 07], supervised feature learning algorithm [Becker 13], and baseline MPP. For a fair comparison, we set the threshold values of the competing algorithms [Talbot 07,Becker 13] to obtain the closest recall scores to the proposed algorithm in Figure 3.13, 3.14 3.15, and 3.16. Blue pixels denote perfectly matching regions as compared with the ground-truth. Green and red pixels show over-detected and under-detected results, respectively.

The main strength of the proposed algorithm is that it ensures stable performances for all datasets without any parameter estimation procedure. The proposed algorithm extracts the most salient line structures in the input image. On the other hand, the proposed algorithm suffers from under-detection when the width of the line structure is varying, for example, see the result for the retina. Such drawback can be overcome if we introduce an additional parameter for width of the line segment in our MPP model.

Conclusions

We propose a novel MPP model to reconstruct curvilinear structures with an optimal set of line segments. For the data likelihood, the proposed density function computes rotated gradient statistics and morphological profiles. SVM reduces the sampling space according to the initially classified areas and defines statistics of the curvilinear features. Prior distribution of the proposed MPP model define interactions of the local configuration in terms of coupling energy states and overlapping areas. We have presented a new optimization scheme which is not biased by the parameter selection in the MPP model. We used an RJMCMC sampler with delayed rejection kernel for different hyperparameter vectors to obtain multiple line hypotheses. The line hypotheses span a feasible sample space, so that the final solution interprets underlying curvilinear structures more faithfully. We have shown line detection results on a wide rage of datasets, and compared the performances of the proposed method with morphological filtering [Talbot 07], supervised learning [Becker 13], and baseline MPP method. The whole optimization process is friendly designed with respect to a parallel implementation; thus, the computational time can be further reduced by applying the parallel Monte Carlo sampler [Verdié 14].

Chapter4

A SSUME that the entire curvilinear structure can be decomposed into many straight line segments. We learn a ranking function that predicts the shape similarity between individual line segments and the latent curvilinear structure. We propose an orientation-aware curvilinear feature descriptor for the learning system, which consists of the image gradient profiles and morphological profiles. Structured Support Vector Machine (SSVM) [Tsochantaridis 05] is employed to obtain weight coefficients for the proposed feature descriptor using the training dataset. To infer the structural information, we build an undirected and weighted graph based on the output rankings of the line segments. We then reconstruct the coarse curvilinear structure by exploring paths of vertices maximizing their distance in the graph. Unlike the previous approaches such as segmentation of the curvilinear structure [Becker 13] and centerline detection [Sironi 14], the proposed algorithm can provide the topological importance level of the curvilinear structure. Figure 4.1 compares results of curvilinear structure segmentation [Becker 13], centerlines detection [Sironi 14], and the proposed method.

To apply a learning system, we need an accurately annotated ground truth for training datasets, where the ground truth is manually labeled by a human expert to classify pixels on the curvilinear structure. It is difficult to secure enough quantity of ground truth images with high quality in that the manual segmentation of curvilinear structure is time consuming work with great efforts. To address this problem, we evaluate a line segment using the shape similarity measure whether it belongs to the latent curvilinear structure. We minimize corruptions of the learning process due to the erroneous annotations.

We also study the graph theory for shape simplification of the complex curvilinear structure. Although the centerline is able to encode the scale of the curvilinear structure, it is insufficient to interpret the topological level of the curvilinear structure. The structured output rankings of the line segment are used to draw the latent curvilinear structure in the order of topological importance. Thus, we create a weighted graph to reconstruct the curvilinear structure with consideration for such topological order. We search for the longest path which minimizes the cumulative rankings of the line segments on the path.

The main contributions described in this chapter are summarized as follows:

• We propose an orientation-aware curvilinear feature descriptor for the ranking learning system;

• We learn a ranking function to predict the correspondence between line segments and curvilinear structure;

• We reconstruct the curvilinear structure with topological order based on graph theory; and

• We test and validate the proposed algorithm on numerous datasets containing micro and macro scales of curvilinear structures in the nature. Das et al. [Das 11] employs SVM to train features and applies skeletonization technique. The overall procedure to detect curvilinear structure is similar to the proposed algorithm. In detail, the proposed algorithm aims to reconstruct curvilinear structures which are latent in different types of image. Thus, we exploit curvilinear features which can encode structural attributes, i.e., image gradient and morphological features. On the other hands, Das et al.'s method targets on a specific application (satellite imaging) so that its feature vectors mainly consist of photometric cues such as RGB histogram features. Use of photometric cue is weak for detecting curvilinear structures with noisy background textures. Also, the proposed progressive path reconstruction algorithm can visualize the topological features in different levels of detail. However, medial axis transform (MAT) employed in [ Das 11] requires an additional process to prune the redundant outcomes which are not associated with the latent curvilinear structure. Unsalan and Sirmacek [Unsalan 12] proposed a road network detection algorithm based on locally predicted probability density function. Thus, the algorithm consists of several refinement steps to obtain the final results, where each step requires manual setting of threshold probability. Unsalan and Sirmacek's method is limited to test many different types of images containing curvilinear structure. On the other hands, the proposed model finds globally optimized line configuration and is capable to detect curvilinear structures on retinal images, biological images, aerial images, and asphalt surfaces.

The rest of the chapter is organized as follows: Section 4.1 reviews the machine learning algorithm on which the proposed method based. Section 4.2 provides the outline of our method. Section 4.3 proposes an orientation aware curvilinear feature descriptor. Section 4.4 explains how we learn the ranking function. Section 4.5 develops a graphical model for inference of curvlinear structure. Section 4.6 shows experimental results on different types of datasets. Finally, Section 4.7 concludes this chapter.

Revisited machine learning systems

We want to extract the structured information of the line networks that appear within homogeneous background textures. Since the information embedded in a single pixel is limited to infer the latent structure, Markov Random Fields (MRF, [Geman 84, Geman 90, Besag 86, Sziranyi 14]) or Conditional Random Fields (CRF, [Lafferty 01, Kumar 06]) based models have been developed to enforce the label consistency on the pre-organized output space. e.g., neighborhood pixels are assumed to have the same label with the high probability when minimizing the cost function. However, the topology of the curvilinear structure composed of line objects is too intricate to be applicable with such design approaches. It is also impossible to specify the topology of the complex line networks with a few parameters. Instead, we employ a machine learning framework, SSVM proposed by [Tsochantaridis 05], to obtain undiscovered distribution which contains the structured relation of input-and-output pairs. The SSVM is often comparable with the graphical models which contain the pairwise struc-ture. Therefore, the SSVM framework is favored to solve computer vision problems, such as object segmentation [Bertelli 11, Lucchi 12, Kim 14], pose estimation algorithm [ Yang 11], and multi-class classification [Mittal 12]. In this work, we use the SSVM framework to train a ranking function that predicts the compatibility with tokenized line segments and the latent curvilinear structures. Finally, the structured output rankings of the line segments are used to estimate the local orientation and to reconstruct the latent curvilinear structures as an ensemble of high ranked line segments. In the remainder of this section, we review Support Vector Machines (SVM) [Vapnik 95] and structured SVM on which the proposed algorithm is based.

Support Vector Machines

The goal of machine learning system is to obtain a function h : X → Z which maps a feature space x ∈ X to an output space z ∈ Z. Using a given training dataset D = {(x i , z i )} K i=1 consisting of input and output pairs, a supervised machine learning algorithm evaluates the quality of a hypothesis ẑ = h(x). Consider a classification problem, we want to find an optimal hypothesis that divides data points into disjoint areas. Figure 4.2 provides a toy example to compare the quality of a hypothesis which is generated by a linear separator denote by black thick line in this figure. Intuitively, we prefer to employ the separator in the bottom row due to the input noise tolerance. Let add a new input point (denoting x in the second columns of the figure) to test and see the behaviors of the two hypotheses. While the input point is spatially closer to the group of red circles than that of blue triangles, the hypothesis of the first row classifies the point as the group of blue triangles. Furthermore, the input point can be perturbed by noise, so that the decision is less reliable at the region near to the boundary. In other words, the noise tolerance can be measured by the distance between a data point and the decision boundary. Thus, the maximum noise tolerance (margin) of the hypothesis is computed using a data point which is the closest to the decision boundary. This is the key idea of the SVM classifier [Vapnik 95, Joachims 99] which maximizes margin for separating given data points. Now, we develop a binary SVM classifier in detail. Suppose that we have K training data points {x 1 , . . . , x K } ⊂ R d and the corresponding labels are z i ∈ {-1, 1}, then a hypothesis is defined as h : R d → {-1, 1}, where d denotes the dimension of the input feature space. Algebraically, a hypothesis is specified as a linearly separable hyperplane w x + b = 0, where w and b are the normal vector and the bias, respectively. So, for classification, we simply check the sign of a given input data point x with respect to the hyperplane:

h(x) = sign(w x + b) (4.1)
Figure 4.2: For the same data points, we compare two linear separators. Intuitively, the separator in the bottom row is safer in that its decision boundary can tolerate a large noise. The maximum noise tolerance of a separator is called its margin corresponding to the thickness of the gray areas in this figure. Thus, SVM is designed to find a hypothesis with the maximum margin for the given training data points.

Let x be a point on the hyperplane. Geometrically, the distance from an arbitrary point x to the hyperplane h is computed by the projection of a vector (xx ) onto w, which is given by

dist(x, h) = |w x -w x | w = |w x + b| w , ( 4.2)
where w x + b = 0. We assume that z i (w x i + b) ≥ 1 to separate all data points with enough margin. It is also called hard margin constraint. Using the fact that z i (w x i + b) > 0, the distance of a training data point x i to h can be written as

dist(x i , h) = z i (w x i + b) w . (4.3)
As we discussed in the beginning of this section, the maximum margin is associated with the distance of the nearest data point to the hyperplane. That is

min i=1,...,K dist(x i , h) = min i=1,...,K z i (w x i + b) w = 1 w . (4.4)
We then formulate the following constrained optimization problem to learn an optimal hyperplane which separates binary classes with the maximum margin:

min b,w 1 2 w w subject to z i (w x i + b) ≥ 1 for i = 1, . . . , K (4.5)
The points which are in the uncertain areas to classify contribute to maximize the margin. In this sense, we conclude that support vectors w are obtained by the nearest points to the hyperplane. The solution of (4.5) is feasible for a duality of a convex optimization problem [ Boyd 04]. More precisely, we solve the following convex optimization problem with Lagrangian multipliers α i 's: While the hard margin constraint is a good strategy for linearly separable dataset, it also induces the over-fitting problem. Moreover, if the given training dataset includes an outlier, decision bound can be dramatically changed (see Figure 4.3). Note that, in many cases, real data are not linearly separable so that the hard margin constraint makes the hypothesis of SVM unpredictable. Thus, we introduce slack variables ξ i ≥ 0 to consider some violated data points which are on the marginal areas (see Figure 4.4). We modify the objective function of SVM (4.5) to take into account the amount of margin violations as the sum of slack variables such as:

L(b, w, α) = 1 2 w w + K i=1 α i (1 -z i (w x i + b)). (4.6)
min b,w,ξ 1 2 w w + C K K i=1 ξ i subject to z i (w x i + b) ≥ 1 -ξ i ; for i = 1, . . . , K (4.7)
where C is a parameter to control the relative importance of total sum of margin violation. We expect that the second term should be minimized to obtain the smallest classification errors. Thus, a large C value induces a hyperplane with thin marginal areas. This modified model is also called soft margin SVM which permits some points to invade the marginal areas.

Structured Support Vector Machines

Despite of the success of SVM models for classification, it is hard to interpret the meaning of the obtained SVM model parameters except that the support vectors are geographically proposed to separate the given training data with the maximum margin. Furthermore, to solve multi-class classification problem with binary SVM classifier, we have to assume the statistical independence between features. More precisely, we isolate one class from the others and then iteratively train SVM classifier to obtain hypothesis of each class. In this scenario, we fails to exploit the correlation information which is latent in the data.

Researchers have developed structure prediction models to deal with complex and interdependent features [Taskar 03, Tsochantaridis 04, Joachims 09]. The joint feature of input-and-output pairs Ψ(x, z) : X ×Z → R N is used to learn structured information. It generates a vector by applying structure z to the given input feature x. For instance, to learn the language structure (grammar), x and z correspond to sample sentence and grammatical structure, respectively. For a ranking learning system, x refers to a list of documents and z is the associated rankings of the documents in the list. A scoring function H : X × Z → R is defined as the inner product of a model parameter w ∈ R N and the joint feature Ψ(x, z), which measures compatibility of the input feature and predicted structure:

H(x, z; w) = w Ψ(x, z). (4.8)
Then, we predict output structure for given model parameter w, which maximizes the scoring function:

ẑ = argmax z∈Z H(x, z; w). (4.9)
For this purpose, we assume that a score with correct structure should be higher than one with incorrect structure. We formulate an optimization problem to predict underlying structure as follows:

min w,ξ 1 2 w w + C K K i=1 ξ i subject to w Ψ(x i , z i ) -w Ψ(x i , ẑi ) ≥ ∆(z i , ẑi ) -ξ i ; ∀ẑ i , ẑi = z i , i = 1, . . . , K, (4.10)
where ∆(z i , ẑi ) denotes a loss function which computes the amount of the loss due to the prediction of ẑi comparing to the ground truth z i . Therefore, for all training data, we constrain that difference between the scores computed by correct structure z i and incorrect structure ẑi is always greater than a required margin. Due to the number of constraints, we need an efficient algorithm to solve the given optimization problem (4.10). Plane-cutting algorithm [Joachims 09] iteratively adds constraints until the prediction error reaches within a pre-defined tolerance value . The algorithm begins to solve (4.10) without constraints, and finds the most violated constraints against for each sample data and the corresponding prediction. For constraints which violate the required margin, plane-cutting algorithm adds ẑ into a working set W i and then solve (4.10) with respect to the updated constraints. This process is also summarized in Algorithm 4.

Algorithm 4 Cutting plane algorithm [Joachims 09]

1: Inputs:

D = {(x i , z i )} K i=1 ∼
a training dataset consisting of input-and-output pairs C ∼ controls a relative importance of regularization and regression terms ∼ a tolerance 2: Output:

w ∼ an optimal model parameter 3: W i ← ∅ for all i = 1, . . . , K 4: while no W i has changed during iteration do 5:

for i = 1, . . . , n do 6: H(x, z; w) = ∆(z i , z) + w Ψ(x i , z) -w Ψ(x i , z i ) 7: Compute ŷ = argmax z∈Z H(x, z; w) 8: Compute ξ i = max{0, max z∈W i H(x, z; w)} 9: if H(x, ẑ; w) > ξ i + then 10: W i ← W i ∪ {ẑ} 11: w ← optimize over W = i W i 12: end if 13:
end for 14: end while

Overview of the proposed algorithm

In this section, we define notations to represent curvilinear structures, and provide an overview of the proposed algorithm. Assume that an image I contains a curvilinear structure. We denote the latent curvilinear structure g : I → {0, 1} for any pixel x:

g(x) = 0 if x is on the curvilinear structure, 1 background. (4.11)
This function is also employed to create a ground truth map, which is manually labeled, for the machine learning system and the performance evaluations. Instead of classifying each pixel whether it belongs to the curvilinear structure, we gather line segments that agree with the latent curvilinear structure of the given image. Hence, the proposed algorithm aims to find a set of line segments s ⊆ S highly corresponding to the curvilinear structure. We define a line segment s ∈ S as a tuple of the center coordinate x, fixed length , fixed thickness τ , and quantized orientations

θ. That is s = (x, , τ, θ) ∈ I × R 2 × Θ, where Θ is set of elements in [0, 2π[.
To train a ranking function with SSVM, we propose a novel curvilinear feature descriptor φ : S → R N that assigns a feature vector φ(s) ∈ R N to each line segment s, where N denotes the dimension of the feature space. Specifically, the feature de-scriptor consists of the morphological profiles φ M to detect the expected length of the structure and the steerable filtering responses φ S to accentuate directional properties (Section 4.3.1). To take into account the directional information of each line segment, we propose an orientation-aware feature vector φ (s) that permutes the elements in the steerable feature vector of φ S (s) according to orientation θ of the given line segment s (Section 4.3.2). Also, we design a spatial neighborhood system for the line segments to improve the consistency of output ranking values (Section 4.3.3). A model parameter w ∈ R N is obtained via machine learning system (SSVM), which determines the relative importance of elements in the curvilinear feature vectors. The output rankings of line segments z ∈ Z are associated to the plausibility for each line segment s to belong to the ground truth map g. For example, a higher ranking is assigned to the given line segment if it fits well to the latent curvilinear structure (Section 4.4). Figure 4.5 summarizes the process how we learn a model parameter of the proposed curvilinear feature descriptor with the given training dataset.

With the rankings of the line segments, we build an undirected and weighted graph G = (V, E) to infer the curvilinear structure with a set of elementary paths on the graph G. Finally, we iteratively reconstruct the curvilinear structure by collecting paths between remote vertices in the graph (Section 4.5).

Orientation-Aware Curvilinear Feature

Feature extraction

This section is devoted to compute the curvilinear feature descriptors that are used for the inputs of the ranking learning system. The evidence of being the curvilinear structure is measured by the homogeneity of the pixel intensities and the directional responses on the path. We extract curvilinear features using the morphological profiles φ M : S → R 4 and the steerable filter responses φ S : S → R 24 . Figure 4.6 depicts an example of the curvilinear feature maps. We refer to Chapter 2 for implementation details of the morphological operation and the steerable filters.

Before applying the filtering operations, we normalize the training and test images to remove the effects of various illumination factors:

Ĩ = 1 1 + e -I-β α , ( 4.12)
where α = max(I) -min(I) represents the range of intensity values, and β = E[I] is the sample mean of the image.

To obtain the morphological profile, we apply the path opening operator [ Talbot 07] with different length parameters L. The path opening operator brightens pixels belonging to the linear shaped objects while it removes other types of structural elements.

As the length parameter of the path opening operator is increased, we obtain a map which emphasizes longer sequence of pixels in the image.

The steerable filtering responses are associated with the directional characteristics of the given line segment. Let F θ be a steerable filter that accentuates orientation θ. The steerable feature map is computed by the convolution of F θ and the normalized image Ī. We construct the baseline steerable filtering feature vector by stacking such steerable filtering responses with increasing order of orientation. That is

φ S (x) = F θ 1 * Ĩ (x) F θ 2 * Ĩ (x) • • • F θ k * Ĩ (x) , (4.13)
where In this work, we use 4 different length parameters to create the morphological profiles and 24 filter banks to build the steerable feature maps. Specifically, the steerable filter banks consist of 3 basis kernels with respect to 8 different orientations. Figure 4.7 shows the steerable filter banks which are employed in this work.

θ 1 < θ 2 < • • • < θ k and k = |Θ| is the number of quantized orientations in Θ ⊂ {0, 2π}.
Finally, the feature vector of the given line segment φ(s) is generated by the concatenation of the feature vectors of φ M (s) and φ S (s). We accumulate the values of the feature maps on the line segment feature and rescale those values according to length of the line segment :

φ(s) = φ M (s) φ S (s) = 1 x∈s φ M (x) φ S (x) (4.14)
Figure 4.8: The baseline steerable feature vector φ S (s i ) is created by steerable filtering responses (F θ * I) as the increasing order of orientation. To take into account the direction of line segment s i , we multiply permutation matrix P θ i by the baseline φ S (s i ).

Permutation of feature maps

As the result of ranking learning system, we obtain a model parameter w that determines the relative importance of each element in the feature vector. If the training set consists of a particular orientation, the corresponding parameter has a larger value than others. However, such situation can induce an over-fitting. To ease this issue, we also exploit the own directional information of the line segments. For this purpose, we permute the order of elements in the baseline steerable feature vector φ S (s) according to the orientation of the given line segment s.

Let P θ i be a permutation matrix that shuffles the order of the feature vector according to the given orientation θ i :

P θ i = e i e (i+n)•k e (i-n)•k • • • e (i+ k /2 )•k , n = 1, 2, • • • , k /2 -1,
where e i denotes a row vector with 1 in the i-th element and 0 in the others. We define a circular modulation operator • which is given by:

i • k = ((i -1) mod k) + 1. (4.15)
For instance, P θ 4 for k = 8 is defined as

P θ 4 =              
e 4 e 5 e 3 e 6 e 2 e 7 e 1 e 8

              =              
0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 To compute the orientation-aware feature vectors φ (s), we apply the permutation matrix P θ to re-order the elements of the baseline feature vector φ(s). For a better understanding, we give a toy example in Figure 4.8 how to shuffle the baseline feature vector for the given orientation. Suppose that the orientation of line segments s 2 and s 5 are associated with θ 2 and θ 5 , respectively. Then, the orientation-aware feature vectors φ (s 2 ) and φ (s 5 ) are given by φ (s 2 ) = φ M (s 2 ) P θ 2 φ S (s 2 ) , φ (s 5 ) = φ M (s 5 ) P θ 5 φ S (s 5 ) .

              . (a) (b) (c) (d)
Now the elements of the orientation-aware feature vectors are aligned in the order of dominant directional response for the given line segments. Therefore, the model parameter w of the learning system properly works to measure compatibility of input line segments and the given image data. Figure 4.8 compares the baseline feature vector with the orientation-aware feature vector for the given line segment s i .

Spatial grouping of the features

To improve segmentation performance, the piecewise smoothness assumption on the image grid is often used [Geman 84, Boykov 04], e.g., 4-connected neighboring pixels have the same label with a high probability. While the thin line structure shows the internal similarity, such prior model also encourages background textures to be classified as the same segment so that it produces discontinuities on the curvilinear structure. Instead, we propose a new topology to consider the neighborhood system of line segments. The proposed neighborhood system is designed to realize a smoothly connected curvilinear structure. Thus, to enhance the spatial coherence of the output rankings, the neighboring set consists of shifted line segments to the extendable direction and shifted-and-rotated line segments with a small curvature. (see Figure 4.9 (a)). We define the spatial grouping operator φ as:

φ(s i ) = j∈N i ω j φ (s j ) j∈N i ω j , ( 4.16)
where N i denotes the neighboring set of the line segment s i and ω j is the corresponding weight factor (Figure 4.9 (a) details the values of ω j used in this work).

Figure 4.9 compares the output rankings obtained by individual feature vector and spatially grouped feature vector. When we learn a ranking function without the spatial feature grouping operation, the output rankings are less reliable due to the corrupted background texture. On the other hand, the output rankings with the spatial feature grouping increase the rankings of line segments corresponding to the latent curvilinear structures.

Learning

In this section, we develop a ranking learning system based on SSVM [ Tsochantaridis 04]. Our goal is to learn a function h : S → Z that predicts structured output rankings of the line segments. For the setup of machine learning system, a training dataset D = {(s i , z i )} K i=1 is required which consists of the inputand-output pairs. Assume that we have a list of line segments s = {s 1 , . . . , s K } ⊆ S and the relevant ranking values z = {z 1 , . . . , z K } ⊆ Z, z i ∈ [0, 1]. The line segments s for training are independently sampled from the images of the training dataset. To determine the corresponding rankings z of the input line segments is a tricky part in this work. Recall that the function g represents a binary map to classify pixels on the latent curvilinear structure. We use the binary map as the ground truth of the curvilinear structure. However, the ground truth g cannot encode the shape information of the line segments (length, orientation, and thickness). In other words, line segments proposed at pixel x could be irrelevant with the true curvilinear structure due to the shape information even though the pixel is on the curvilinear structure,i.e., g(x) = 0. We evaluate the shape dissimilarity between the given line segments s i and the corresponding patch of the ground truth map g i :

z i = 1 -score(s i , g i ).
(4.17 In this work, we employ F 1 measure [Martin 04, Arbelaez 10] which is defined as

score(s i , g i ) = F 1 (s i , g i ) = 2P i R i P i + R i , ( 4.18)
where P i and R i denote precision and recall, respectively [Philbin 07, Manning 08], which are associated with the line segment s i . Figure 4.10 shows an example that computes the shape similarity score of a line segment and the corresponding ground truth patch. If s i lies on the true curvilinear structure g i , the corresponding ranking value z i is closer to 0; otherwise s i is on a non-curvilinear structure, z i is set to 1.

We use a joint feature map of input-and-output pairs Ψ(s, z) : S × Z → R N to obtain structured output rankings. A scoring function h : S × Z → R is defined as the inner product of the joint feature vector and a model parameter w. For a given line segment, a prediction ẑ = h(s; w) is performed by finding rankings that maximize the scoring function: We define the joint feature map Ψ(s, z) to obtain the decreasing order of rankings among input line segments:

Ψ(s, z) = ij z ij φ(s i ) -φ(s j ) , ( 4.20)
where the ranking matrix z ij indicates the pairwise order of the objects s i and s j in terms of shape dissimilarity measure (4.17). In other words,

z ij = +1 if z i > z j , -1 otherwise.
The parameter of the learning system w ∈ R N is obtained by minimizing the following constrained objective function:

min w,ξ 1 2 w w + C K K i=1 ξ i subject to w Ψ(s, z) ≥ w Ψ(s, ẑ) + ∆(z, ẑ) -ξ; ∀ẑ, ẑ = z (4.21)
where the first term of (4.21) penalizes the complexity of the solution (regularization) and the second term measures the relaxed loss of the system to fit the parameter into training data (regression). C controls the relative importance between the two terms. ∆(z, ẑ) is a loss function that counts the number of incorrectly ranked pairs. To solve this optimization problem, we use the cutting plane algorithm ([Joachims 09], see Algorithm 4 in Section 4.1.2).

We have obtained a ranking function h(s; w) that arranges the line segments in the order of correspondence with the underlying curvilinear structure. Specifically, the linear combination of the model parameter w and the proposed feature descriptor φ(s), ẑ = w φ(s), is used to evaluate the line segment s. At each pixel grid, we create k line segment candidates, where k denotes the number of orientations. We then estimate the local orientation using the ranking function:

î = argmax i∈{1,...,k} w φ(s(θ ← θ i )), then θ max = θ î, (4.22)
where s(θ ← θ i ) denotes a candidate of the line segment which updates its orientation for θ i whereas the rest of the parameters (center locations, length, and thickness) are fixed. Moreover, the ranking function is employed to indicate a pixel on the curvilinear structures. For this purpose, the output rankings are computed at each pixel grid using the line segment maximizing local orientation responses via (4.22).

Inference

In this section, we aim to develop a graphical model to infer the topology of the curvilinear structure. The main motivation to represent the curvilinear structure as a graph is to use an optimal number of pixels to reconstruct the curvilinear structure without redundancies. Also, we want to see the topological features of the structure at different levels of detail.

For the sake of convenience, we recall several terminologies related to graph theory before developing the proposed algorithm. We consider the graph G = (V, E) where V is the set of pixels. Each pixel being given its image coordinates, and pixels u = (i, j) and v = (m, n) are adjacent if and only if |i -m| ≤ 1 and |j -n| ≤ 1. Moreover, we assign a weight for each edge {u, v} ∈ E as w(u, v) = uv ẑu+ẑv 2 , where uv denotes the Euclidean distance between pixel coordinates of u and v. ŷu and ŷv are the output rankings at u and v, respectively. A path P in the graph G is a sequence of distinct vertices such that consecutive vertices are adjacent. The length of P is the sum of the weight of its edges, and the distance dist(u, v) between two vertices u and v is the minimum length of a path from u to v. The eccentricity ecc(v) denotes the maximum distance from the node v to a node u ∈ V , i.e., ecc(v) = max u∈V dist(u, v). The diameter diam(G) of G equals max v∈V ecc(v), i.e., it is the maximum distance between two nodes in G.

Intuitively, to infer the topology of the curvilinear structure, we look for long shortest paths in the subgraph G of G induced by the nodes with high rankings. More precisely, our algorithm computes a diameter of G , i.e., a shortest path P with length diam(G ). This path P is added in the curvilinear structure inference, then the path P is contracted into a single node (equivalently, the weight of all edges of P become 0), and the process is repeated till the diameter of the subgraph is smaller than pre-defined path length min . The entire procedure is summarized in Algorithm 5.

P ← P ∪ Q 12:
w(u, v) ← 0, ∀u, v ∈ Q, All edge weights on the path Q become 0 13: end while

The first part of the algorithm is to extract the subgraph G from the whole set of pixels. This subgraph should contain the curvilinear structure and be relatively small for an efficient computation. For this purpose, from the training datasets, we first compute the average proportion ρ of pixels being part of the curvilinear structure. Then, the binary segmentation map is induced by ρ|I| pixels according to the output rankings, related to a high ranking. However, due to the noises and the choice of a small ρ, the binary segmentation map may be disconnected so that some parts do not reflect the connected components of the curvilinear structure. Moreover, the binary segmentation map only encodes the center coordinates of the line segments (see Figure 4.11) while the output rankings contain the shape information of the line segments. To overcome these issues, we disjoint the binary segmentation map into line segments, and compute the dissimilarity scores (4.17) between these line segments and the corresponding patches of the binary segmentation map. We fill up areas corresponding to the line segments with those of dissimilarity scores. Specifically, for overlapping areas between adjacent line segments, we assign the average value of theirs (see Figure 4.12). Thus, the subgraph G is built upon a subset of pixels V ⊂ V which corresponds to the dissimilarity score map. We also update the weights for each edge {u, v} ∈ E in the subgraph as w(u, v) = uv bu ẑu+bv ẑv 2 , where b u and b v denote the dissimilarity scores at u and v, respectively.

The time-consuming part of the proposed algorithm is the computation of a diameter of G . Rather than computing all pair distances (which requires a linear number of application of Dijkstra's algorithm), we use an efficient heuristic algorithm called 2-sweep algorithm proposed by [Corneil 01]. The 2-sweep algorithm randomly picks a vertex t in G , then performs Dijkstra's algorithm from t to find a node u at the maximum distance from t, i.e., dist(u, t) = ecc(t). Then, it computes (using Dijkstra's algorithm) a path from u to a node v at the maximum distance from u. The length of the second path (from u to v) is a good estimation of diam(G ). Note that this algorithm is able to compute the exact diameter if the graph has tree structure [Borassi 15]. Topologically speaking, most of latent curvilinear structures are very close to trees, so that the use of the 2-sweep algorithm is well adapted. For a better understanding, we schematically explain the intermediate steps of the proposed curvilinear structure inference algorithm in Figure 4.13.

Experimental results

In this section, we first discuss the parameters of the proposed algorithm according to the datasets. We then compare the quantitative and qualitative results of the proposed algorithm and those of competing models proposed by [Frangi 98], [ Law 08], [Becker 13], [Sironi 14], and [Jeong 15a].

Parameters

The proposed algorithm requires few parameters to compute the curvilinear feature descriptor φ for defining physics of the line segment,i.e., length and thickness τ , and to train a ranking function via SSVM. For the morphological profiles φ M , we apply the path opening operator 4 times with different length parameters L. We linearly increase such length parameters by adding η from the minimum length L * . Thus, L ∈ {L * + nη|n = 0, 1, 2, 3}. In the case that the curvilinear structure appears darker than surrounding pixels, we invert the input images before applying the morphological operation. This is because the path opening operator assumes pixels on the curvilinear path to be brighter than neighboring pixels. For the steerable features φ S , we perform convolution of the steerable filters with the normalized input images. The steerable filter banks consist of 3 basis kernels with k = 8 different orientations, where the size of the filters is 21 × 21 (see Figure 4.7). To build the basis kernels, we superpose the elongated Gaussian kernels using different aperture and smoothing factors. The aperture value determines the space between the Gaussian kernels with respect to the minor axis. We control the aperture size of the Gaussian kernels as a = {6, 3, 1} to take into account the thickness of the curvilinear structure. The smoothing factor of the outer kernels σ 2 o is set to 2.5 and those of inner kernels σ 2 i is set to 1.5. The elongation factor γ is fixed as 2.5 for all experiments. Note that the steerable filters are symmetric with respect to the major axis. Also, we assume that the length of line segment should be larger than its thickness, > τ , and these two values depend on the dataset. For SSVM training, we sample 500 positive and negative query line segments, respectively, for every training image. C which controls the relative importance of regularization and regression terms in (4.21) is set to 0.01 for all datasets. We choose the set of parameters for the SSVM training via 3-fold cross validation [Kohavi 95] which maximizes the average F 1 score of the training set.

We test our curvilinear structure model on the following public datasets:

• DRIVE [Staal 04]: This 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, re-spectively. We build the morphological profiles with (L * , η) = (40, 10) and use ( , τ ) = (5, 3).

• RecA [Jeong 15b]: We collect electron microscopic images of RecA proteins on DNA which contain filament structure. We use 4 training images and 2 test images. We set (L * , η) = (40, 10), and ( , τ ) = (5, 3), respectively.

• Aerial [Sironi 14]: The dataset contains 14 remote sensing images of road networks. We select 7 images for the training phase and 7 images for the testing phase, respectively. (L * , η) is set to (50, 50) and each line segment consists of = 15 and τ = 4.

• Cracks [Chambon 10]: Images of the dataset correspond to road cracks on asphalt surfaces. We use 6 images to train and test the algorithms on different 6 images. Parameters for the morphological profiles (L * , η) are set to (10, 5), and for the physics of line segments ( , τ ) are set to (5, 3).

To create subgraph G for the curvilinear structure inference, we compute the average proportion ρ of pixels on the curvilinear structure by observing the training images. We obtain ρ values for DRIVE, RecA, Aerial and Cracks dataset as 8.63%, 4.99%, 8.16%, and 4.3%, respectively. The minimum path length ˆ is adaptively set to 40, 30, 80, and 30 for DRIVE, RecA, Aerial and Cracks datasets, respectively, to obtain the best performance.

Evaluations

If we use pre-computed feature descriptors, sorting the given line segments takes less than a seconds. The proposed algorithm progressively reconstructs the curvilinear structure by adding a long path on the subgraph. Figure 4.14 and Figure 4.15 show the intermediate steps of the proposed graph inference algorithm for Drive and Aerial dataset, respectively. Unlike the previous models, the proposed algorithm is able to show different levels of detail for the latent curvilinear structure. Such information, which visualizes shape complexity of the curvilinear structure, cannot be retrieved by tuning the threshold. In practice, a few number of iterations is required to converge the proposed 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. Computing curvilinear feature descriptors is the most time consuming part. To compute all possible line segments in 800×500 pixels image, it takes 5 minutes -we expect that parallel implementation can dramatically reduce the computation time. To train SSVM, we use 1000 examples and it takes a minute. We visually compare the performance of the proposed algorithm with the competing algorithms. Figure 4.16, 4.17, 4.18, and 4.19 show the results from DRIVE, RecA, Aerial, and Cracks dataset, respectively. Although the proposed algorithm is the most suitable to show the topological importance level of the curvilinear structures, its intermediate result (dissimilarity score map) is comparable to the results of curvilinear segmentation algorithms.

For the quantitative evaluation, we provide precision-and-recall curve of the proposed algorithm (dissimilarity score map) and the state-of-the-art models in 4.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. To provide the average performance of each dataset, we compute the interpolated average precision values at the fixed recall values, where the method is detailed in [ Manning 08]. However, the evaluation with precision and recall is unable to measure the structural features. We also provide the average proportion of pixels to represent the curvilinear structures and the corresponding F 1 scores in Figure 4.21 and Table 4.2. The graph shows that the proposed algorithm efficiently draws the curvilinear structures using a smaller number of pixels than the other algorithms. In Table 4.2, we observe that the proposed algorithm outperforms than competing algorithms for all datasets, except Drive dataset. It is because the proposed algorithm uses the fixed thickness τ to describe line segments whereas the images of Drive dataset consists of blood vessels with varied thicknesses. The performance could be improved if we consider the thickness of the line segment as a learnable variable.

It is worth to mention the over-fitting problem in curvilinear segmentation task

Conclusions

This chapter proposes a new curvilinear structure reconstruction algorithm based on the ranking learning system and graph theory. The output rankings of the line segments correspond to the plausibility of the latent curvilinear structure. Using an optimal number of pixels, the proposed algorithm provides different levels of detail during reconstruction of the curvilinear structure. More precisely, we learn a ranking function based on SSVM with the proposed orientation-aware curvilinear feature descriptor. The weighted average of curvilinear feature groups improved the consistency of the output rankings. In this chapter, we define a novel graphical model that infers the curvilinear structure according to the topological importance. The proposed algorithm looks for remote vertices on the subgraph which is induced from the output rankings. Across the various types of datasets, our model shows good performances to reconstruct the latent curvilinear structure with a smaller number of pixels compared to the stateof-the-art algorithms. Furthermore, to our best knowledge, this is the first work to provide the topological importance of the curvilinear structure.

I N this dissertation, we proposed curvilinear structure reconstruction models which can interpret various types of natural images. This research topic let us ask a question. "How can we well represent curvilinear structure in general situations?" Geometrically speaking, curvilinear structures exhibit thin, elongated, and symmetric for a particular direction. Also, pixels on the curvilinear structure show slightly different intensity values compared to its surroundings. These are generally accepted assumptions in order to describe curvilinear structure in computer vision and image processing. Thus, oriented image gradient information and morphological filtering responses have been employed to extract such curvilinear features. In this thesis, we specifically assumed that the entire curvilinear structure can be separated into multiple line segments. We look for a set of line segments corresponding to the above mentioned characteristics of curvilinear structure. Each line segment is parameterized as its center pixel coordinate, length, and orientation. We further studied recently developed computer vision techniques, e.g., stochastic modeling, structured information prediction, and graph theory.

In Chapter 2, we first reviewed the existing curvilinear structure reconstruction models to analyze blood vessels, biological images, facial wrinkles, road networks, and defects in the asphalt. Most of the previous works have been designed for specific applications due to the lack of evidence to detect curvilinear structure using local image features and complex topological behaviors depending on the applications. We found a common description of the curvilinear structure based on its geophysical characteristics, such that the curvilinear structure exhibits thin, elongated, and symmetrical structure on images. We obtained corresponding image features via rotated gradient information and morphological image filtering. We also realized that geometric priors are essential to improve the performances by inferring missing information.

We implemented two different curvilinear structure reconstruction models with stochastic model in Chapter 3 and machine learning in Chapter 4. Chapter 3 began to propose a probability density of line segments for given image data. To maximize the proposed density, we employed RJMCMC sampler with delayed rejection. The proposed stochastic model requires a small number of hyperparameters to decide relevant 101 importance between geometric constraints. Instead of estimating the exact values of these hyperparameters, we simulated several Markov chains to obtain line hypotheses. An aggregation algorithm of these line hypotheses is proposed to reduce sampling space and to find the most promising parameter set.

Chapter 4 proposed a new curvilinear structure representation method based on graph theory which exploits output ranking of the line segments. The proposed curvilinear structure inference algorithm visualized topological importance according to the length of the structure. We learned a ranking function to predict the correspondence between line segment and the latent curvilinear structure. The proposed algorithm efficiently reconstructed curvilinear structure with the minimum number of pixels. For each chapter, we compared the quantitative and qualitative performances of the proposed algorithm and of the competing methods.

The proposed algorithms aim to find an optimal set of line segments to represent the latent curvilinear structure. The same curvilinear features (steerable gradient profiles and morphological profiles) are employed in the proposed models. In Chapter 3, stochastic model for curvilinear structure is proposed to localize a set of line segments onto the latent curvilinear structure. It is useful for application to exploit information regarding distribution of the line segments. For instance, to evaluate facial wrinkle condition, the density of the line segments can be a criterion. In Chapter 4, we propose a curvilinear structure reconstruction algorithm based on ranking learning system and graph traversal algorithm. The later proposed algorithm has the merit of preserving topological feature of the curvilinear structure. Thus, it shows robust reconstruction results without discontinuities, where the input image consists of noisy background textures.

Perspectives

Our stochastic model proposed in Chapter 3 often showed disconnected line segments at junctions (see Figure B.2). It is because the prior model is designed to reject large curvatures in the configuration. We can remedy this drawback by employing another layer to represent the curvilinear structure in a different form e.g., Bezier curve networks [Bama 15] or a sequence of pixels. We define a transform function, which is invertible, from a group of line segments to a desired representation. We assume a strong connectivity of the pixels belonging to the same group. To implement this model, we shall design a new transition kernel for MCMC sampler which can jump between different representation layers in a reversible manner.

For our second model regarding machine learning, we can improve the feature descriptor to take into account the symmetric shape. This term will contribute to high- light output rankings at the center of the linear structure. The idea is inspired by the work of [Sironi 14] who exploits distance transform to identify centerlines. Moreover, we consider training non-line shape templates such as bifurcation patterns [Azzopardi 11] in order to localize junctions.

In this thesis, each line segment is parameterized as its center pixel coordinate, length, and orientation. Missing the parameter to describe thickness, the proposed algorithms show performance drops for detection of width-varying curvilinear structures i.e., blood vessels. To overcome this issue, we can add a thickness parameter to describe line segment. We also consider extending the proposed models with hierarchical grouping.

Certainly, the idea involving machine learning technique (employed in Chapter 4) with MCMC optimization (proposed in Chapter 3) is very interesting. We expect that the hyperparameters can be learned and updated at the each iteration of MCMC optimization (coupled with simulated annealing). Specifically, Wick et al. [Wick 11] proposed the ranking learning algorithm called SampleRank which can be combined with MCMC optimization. However, our primitive attempt to apply this method into the proposed cost function has failed due to the following reasons: First, the cost function of the stochastic model is non-convex because of the overlapping penalty term (3.16) that prevents congestion of the line segments in the local configuration. Second, since the groundtruth map is provided by a binary segmentation map, we do not have

A.5 Crack: Defects in the asphalt

Images regarding road cracks contain thin and elongated structures in the asphalt surfaces (see Figure A.5). The dataset provides the real and synthetic images of various resolutions. We select 12 real images in the dataset, which have the same resolution (768 × 512). We use 6 images to train, and different images 6 to test the algorithms. Please note that images are under copyright of the IFSTTAR institute, France. We obtain this experimental data courtesy of the authors [Chambon 10]. donnée afin d'obtenir les directions principales de la structure locale [ Frangi 98]. Le flux optimal orienté (Optimally Oriented Flux en anglais, OOF) mesure la quantité de flux de gradient sortant pour trouver les structures curvilignes [ Law 08]. L'opérateur morphologique recueille des pixels en fonction de la similitude structurelle sur le chemin allongé [ Talbot 07]. Toutefois, en raison de l'absence d'interprétation de forme, ces méthodes fondées sur les caractéristiques locales d'image sont insuffisantes pour reconstruire la structure curviligne sous-jacente. Par ailleurs, des modèles graphiques tels que [González 10,Türetken 13b] définissent des contraintes géométriques dans une configuration locale en minimisant globalement leurs fonctions de coût. Plus précisément, les algorithmes fondés sur des graphes sont initialisés par des points correspondant à des structures curvilignes latentes, puis définissent un trajet qui relie ces points avec un a priori géométrique pour fournir des formes plausibles. Les propriétés géométriques de réseau de lignes sont impliqués comme termes de contraintes quand un problème d'optimisation de l'énergie est formulé. Des modèles stochastiques [Lacoste 05, Jeong 15a] reconstruissent les structures curvilignes en échantillonnant plusieurs segments de ligne afin de maximiser une probabilité a posteriori de données d'image fournie. Des concepts allant de la similarité à la représentation graphique sur des a priori géométriques sont considérés pour définir la connectivité et la courbure des segments de ligne. Récemment, des algorithmes d'apprentissage automatique ont été proposés pour détecter la structure curviligne. Becker et al. [Becker 13] ont appliqué un algorithme pour obtenir un ensemble optimal de bancs de filtres de convolution. Sironi et al. [Sironi 14] ont développé un modèle de régression pour estimer l'échelle (largeur) des structures curvilignes et localiser leurs axes.

Bien que les algorithmes de [Tu 06, Arbeáez 11] détectent des courbes et des lignes, l'objectif est tout à fait différent des techniques de reconstruction des structures curvilignes. Les algorithmes de contour mentionnés ci-dessus cherchent des lignes de contour fermé pour diviser une image en régions significatives. Par conséquent, des fonctions de coût exploitent les signaux globaux de texture en ce que les contours sont associés à des bords saillants près des limites de l'objet. A l'opposé, nous recherchons plusieurs structures curvilignes, qui ne sont pas nécessairement fermées et sont latentes dans une texture homogène. En comparaison avec le contour, les structures curvilignes sont estimées via de subtiles caractéristiques d'images locales. Les similitudes internes de la structure et la conception précise de la forme sont essentiels pour résoudre notre problème.

B.2 Contributions

Dans cette thèse, nous avons comme objectif ambitieux de proposer un cadre général qui peut analyser différents types de structures curvilignes variées. Nous supposons que la structure curviligne peut être décomposée en un ensemble de segments de ligne. Nous cherchons donc un ensemble optimal de segments de ligne qui correspond à la structure curviligne latente. Nous avons proposé deux modèles de reconstruction de la structure curviligne fondés respectivement sur un modèle stochastique avec une technique d'échantillonnage (Chapitre 3) et sur un classifieur à partir d'un système d'apprentissage avec une représentation graphique (Chapitre 4). Pour le modèle stochastique, nous proposons une probabilité a posteriori de segments de ligne pour l'image donnée, et nous optimisons la répartition proposée par un echautillonnage de type Markov chain Monte Carlo (MCMC). Notre travail contribue spécialement à réduire la dépendance pour les hyperparamètres du processus stochastique. À cette fin, nous simulons simultanément plusieurs chaînes de Markov avec différents paramètres et des hypothèses multiples de lignes agrégées. En outre, nous apprenons une fonction de classement qui mesure la comparabilité des segments de ligne donnés et les structures curvilignes latentes. Une méthode de représentation graphique utilisant les classements de sortie des segments de ligne est proposée pour montrer l'importance de la structure topologique curviligne sous-jacente. Nous testons différents types d'ensembles de données, et fournissons des résultats expérimentaux dans chaque chapitre.

B.2.1 Chapitre 2 : Modélisation de la structure curviligne : présentation générale

Dans le chapitre 2, nous fournissons des descripteurs de structure curviligne fondés sur le gradient de l'image et des profils morphologiques. Nous percevons la structure curviligne latente sur la base de l'incohérence des textures de fond et des caractéristiques géométriques. En d'autres termes, une séquence de pixels correspondant à la structure curviligne a des valeurs d'intensité différentes comparé à son environnement, et exhibe une forme fine et allongée. Le gradient de l'image est calculé par rapport à diverses orientations pour détecter des fonctionnalités orientées localement. En outre, nous employons des profils morphologiques pour trouver des pixels représentant les differéntes textures du fond sur un chemin allongé. Dans ce chapitre, nous examinons aussi les travaux antérieurs sur la structure des modèles de reconstruction curviligne tels que la segmentation, la représentation à base de graphes, la modélisation stochastique, et l'apprentissage automatique (machine learning en anglais).

B.2.2 Chapitre 3 : Modèle stochastique

Dans le chapitre 3, nous étudions un modèle stochastique pour la reconstruction de la structure curviligne. Nous proposons une densité de probabilité de segments de ligne pour une image donnée, où cette densité de probabilité est constituée de la vraisemblance et de l'a priori. Pour la probabilité au terme d'attache aux données (vraisemblance), nous exploitons le gradient de l'image pivotée et les réponses des filtres morphologiques pour estimer la présence d'une structure curviligne. Nous approchons les statistiques de caractéristiques curvilignes en utilisant une distribution gaussienne avec un SVM (Support Vector Machine). Nous limitons aussi des comportements géométriques locaux de segments de ligne. Nous encourageons la bonne connexion de segments de ligne, alors que nous rejetons l'encombrement de segments de ligne. Pour optimiser la densité proposée, nous employons un algorithme de type MCMC à saut réversible (RJMCMC) avec un rejet retardé. La principale contribution de ce chapitre est d'éliminer la dépendance de la performance de la sélection des hyperparamètres en utilisant l'optimisation stochastique. Plusieurs chaînes de Markov sont simulées avec différents vecteurs d'hyperparamètres afin d'obtenir différentes hypothèses de ligne. Puis nous agrégeons ces hypothèses afin de réduire l'espace d'échantillonnage. Nous fournissons des expériences pour comparer les résultats de reconstruction de la structure curviligne obtenue par l'algorithme proposée et des algorithmes concurrents de l'état de l'Art.

B.2.3 Chapitre 4 : Inférence de la structure curviligne

Le chapitre 4 traite d'un modèle d'inférence de structure curviligne qui quantifie un score en fonction de l'entrée et de la sortie des informations de paires de segments de ligne. Pour la classification, les méthodes antérieures définissent globalement un seuil pour supprimer des pixels qui montrent une faible probabilité d'appartenir à la structure curviligne. Au cours de ce processus, la plupart des travaux précédents ignorent implicitement les informations corrélées des pixels situés sur la structure curviligne. Nous visons à apprendre une fonction de classement des segments de ligne pour des structures curvilignes données. Le classement de sortie est en mesure d'organiser les segments de ligne dans l'ordre de correspondance avec les structures curvilignes latentes. Plus précisément, un SVM Structuré (SSVM) est utilisé pour déduire les classements de sortie. Nous visualisons ensuite les ordres topologiques de la structure curviligne. Nous étudions la théorie des graphes pour reconstruire le trajet curviligne qui relie les plus reculés des pixels sur la structure curviligne latente. L'algorithme proposé ajoute progressivement les branches les plus fines dans la structure principale jusqu'à ce qu'il trouve la longueur minimale de la structure. Nous montrons des résultats expérimentaux sur divers ensembles de données pour illustrer l'efficacité de l'algorithme proposé.

B.3 Conclusions

Dans cette thèse, nous avons proposé des modèles de reconstruction de structures curvilignes qui peuvent interpréter différents types d'images naturelles. Posons une question relative à notre sujet de recherche : "Comment pouvons-nous bien représenter la structure curviligne dans des situations générales?" Pour répondre à cette question, nous supposons que la structure curviligne peut être représentée comme une collection de nombreux segments de ligne. Nous avons étudié davantage des techniques récemment développées en vision pas ordinateur, par exemple, la modélisation stochastique, la prédiction de l'information structurée, et la théorie des graphes.

Dans le chapitre 2, nous avons d'abord examiné les modèles de reconstruction de structures curvilignes existants pour analyser des vaisseaux sanguins, des images biologiques, des rides du visage, des réseaux routiers et des défauts dans l'asphalte. La plupart des travaux précédents ont été conçus pour des applications spécifiques en raison de l'absence de preuve pour détecter des structures linéaires curvilignes utilisant des caractéristiques de l'image et des comportements locaux topologiques complexes en fonction des applications. Nous avons proposé une description commune de la structure curviligne basée sur des caractéristiques géométriques, telles que la finesse, l'allongement et symétrie sur les images. Nous avons obtenu des caractéristiques de symétrique. Ce terme contribuerait à mettre en évidence les classements de sortie au centre de la structure linéaire. L'idée est inspirée par le travail de [Sironi 14] qui exploite la transformée de distance pour identifier les axes. Aussi, nous pourrions considérer des structures non-curvilignes telles que des bifurcation en V ou en T par exemple, afin de localiser les discontinuités [ Azzopardi 11].

Nous pourrions aussi accélérer les algorithmes proposés en employant la programmation parallèle. Plus précisément, le processus d'optimisation de notre modèle stochastique est parallelisable. Le temps de calcul peut être encore réduit par l'implémentation parallèle de l'échantillonneur RJMCMC [Verdié 14].

xiii 2. 6

 6 The graph-based models reconstruct curvilinear structures by connecting seed points (denoted by green points), which are associated with local maxima of tubularity scores. (a) The log-likelihood of the image based termlog P (I, G|t) is minimized to recontruct path, which shows topological errors (zig-zags and surplus branch). (b) Geometric priorlog P (t) is incorperated to correct such topological erors. (c) Additional relaxation term is required to prevent counting a seed point twice for path reconstruction. Images courtesy of[Türetken 11]. . . . . 2.7 Stochastic models constain the spatial distribution of the line segments, where each line segment s i is represented as a tuple of its center coordinate x i , length i , and orientation θ i . (a) Pair line segments are assumed to be connected if the end-to-end distance d ij of them is relatively small (< ). (b) The angular difference θ ij of two line segments is related to the curvature. A small curvature is preferred in order to reconstruct smoothly connected curvilinear structures. (c) Congestion of line segments will be penalized. . . . . . . . . . . . . . . . . . . . . 2.8 Intermediate results of the learning based SSDT algorithm: (a) input image, (b) centerline, (c) distance transform, (d) a function to be used for learning, (e) output function of Gradient Boosting, and (f) applying NMS on the output function to localize centerlines. Images courtesy of [Sironi 14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 A spatial point process is a mathematical model to analyze spatial point patterns using a finite number of parameters. More precisely, the spatial point process is a random variable which maps an unordered set of points in a finite space. We can compare different point patterns according to the randomness in the space. . . . . . . . . . . . . . . . . . . . . . . . 3.2 Comparison with the realization of point patterns obtain by (a) homogeneous Poisson process, (b) inhomogeneous Poisson process, and (c) repulsive Markov point process. For a finite space, homogeneous Poisson process employs a constant intensity measure, whereas inhomogeneous Poisson process employs intensity measure as a function. Both Poisson processes ignore correlation of points in the space. On the other hand, Markov point process constrains the local interaction of neighboring points using a non-negative density function. Here, repulsive constraint prevents overlapping between points. . . . . . . . . . . . . . 3.3 Line segments (s m , s n ) interact if they are close enough. . . . . . . . . .3.4 For the input image (a), we compare the manually labeled result (b), outcome of the binary SVM classifier (c) and the posterior probabilities (d). The proposed stochastic model (e) refines the weak points of the SVM classifier by sampling line segments with geometry prior constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 For coupling feature c mn We measure (a) the distance of adjacent line segments d mn and (b) the angle differnece between line segments θ mn .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 The MCMC sampler perturbs the current configuration s via subtransition kernels: (a) birth kernel, (b) death kernel, and (c) lineartransform kernel. Schematic diagrams illustrate how the sub-transition kernels change topology of the corresponding configuration, and a state of the Markov chain is gradually updated in order to maximize the proposed probability density function. If a new configuration disagrees with the probability density function, the chain will remain at the previous state as shown in (d). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Given configuration (a), if a line segment s proposed by LT kernel is rejected (b), the delayed rejection kernel searches for the nearest extremes in the rest of line segments (c). An alternative line segment s, which enforces connectivity, will be proposed by interpolation of the retrieved points (d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Given the input image (a), we compute the gradient magnitude (b).

4. 3

 3 Even a single outlier (see black dot) can dramatically change the characteristic of the hyperplane. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Curvilinear structures are observed in different topological features under the various background textures. Left to right images show the latent curvilinear structure in blood vessels, DNA filaments, facial wrinkles, and road cracks.

  (a) Input, I (b) F θ1 * I, θ 1 = 0 • (c) F θ2 * I, θ 2 = 22.5 • (d) F θ3 * I, θ 3 = 45 • (e) F θ4 * I, θ 4 = 67.5 • (f) F θ5 * I, θ 5 = 90 • (g) F θ6 * I, θ 6 = 112.5 • (h) F θ7 * I, θ 7 = 135 •

Figure 2 . 1 :

 21 Figure 2.1: The local orientation of the given image I is analyzed by performing image convolution with rotatable filter banks. Filtering responses accentuate pixels according to the specified orientation. The associated convolution filters F θ are depicted in left bottom of each image.

Figure 2

 2 Figure 2.2: A candidate kernel is generated by the linear combination of shifted anisotropic Gaussian kernels. This approach is employed to analyze the size of underlying curvilinear structure intuitively.

Figure 2 . 3 :

 23 Figure 2.3: We compare convolution filter banks which are built constructively (top row) and analytically (bottom row). Increasing smoothing factor σ is related to estimating width of the structure, it is more likely to remove noises. On the other hand, controlling the aperture value a of constructive filters provides consistent result that highlight pixels on the curvilinear structure according to the width. The associated convolution filters are depicted in left bottom of each image.

Figure 2 . 4 :

 24 Figure 2.4: Toy example to show how the path opening operator finds valid paths on the adjacency graph: (a) an adcaency graph composed with south-north cones, (b) greed nodes are associated with pixels after image binarization, (c) seeking paths, denoted red color, on the graph, and (d) a path is valid if its length is larger than pre-defined length L.

Figure 2 . 5 :

 25 Figure 2.5: Morphological profiles: Path opening [Talbot 07] results with different length parameter L. As the length parameter is increased, the path opening operator removes small objects which are not on a valid path.

1.

  Curvilinear structure segmentation [Frangi 98, Law 08] 2. Graph-based models [González 10, Türetken 13b] 3. Stochastic image models [Lacoste 05] 4. Machine learning based approaches [Becker 13, Sironi 14]

Figure 2 . 6 :

 26 Figure 2.6: The graph-based models reconstruct curvilinear structures by connecting seed points (denoted by green points), which are associated with local maxima of tubularity scores. (a) The log-likelihood of the image based termlog P (I, G|t) is minimized to recontruct path, which shows topological errors (zig-zags and surplus branch). (b) Geometric priorlog P (t) is incorperated to correct such topological erors. (c) Additional relaxation term is required to prevent counting a seed point twice for path reconstruction. Images courtesy of [Türetken 11].

  where Θ is set of elements in [0, 2π[. A probability density function of line segments p(s|I) is proposed to fit line segments into the given image data and to constrain spatial relationship of the line segments. Thus, the probability density function consists of a data dependent term p(I|s) and a prior energy p(s) given by p(s|I) ∝ p(I|s)p(s), (2.19) where s = {s 1 , . . . , s n } denotes a configuration of the line segments. More specifically, to localize line segments onto the curvilinear structure, the data term p(I|s) exploits curvilinear features such as local statistics of pixel intensities [Lacoste 05], Gabor filtering responses [Tu 06], and image gradients [Batool 12, Jeong 14]. The prior energy p(s) defines local interaction of the line segments in terms of the connectivity and the curvature (see Figure 2.7).

Figure 2 . 7 :

 27 Figure 2.7: Stochastic models constain the spatial distribution of the line segments, where each line segment s i is represented as a tuple of its center coordinate x i , length i , and orientation θ i . (a) Pair line segments are assumed to be connected if the endto-end distance d ij of them is relatively small (< ). (b) The angular difference θ ij of two line segments is related to the curvature. A small curvature is preferred in order to reconstruct smoothly connected curvilinear structures. (c) Congestion of line segments will be penalized.

  Figure 2.8 provides the intermediate steps of SSDT algorithm proposed in [Sironi 14].

  Figure 2.8: Intermediate results of the learning based SSDT algorithm: (a) input image, (b) centerline, (c) distance transform, (d) a function to be used for learning, (e) output function of Gradient Boosting, and (f) applying NMS on the output function to localize centerlines. Images courtesy of [Sironi 14].

  Contents 3.1 Revisited Marked Point Process Framework . . . . . . . . 3.2 Stochastic modeling of curvilinear structure . . . . . . . . 3.2.1 Data likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Prior Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Markov chain Monte Carlo sampler . . . . . . . . . . . . . . . 3.3.2 Reversible Jump Markov chain Monte Carlo sampler . . . . . 3.3.3 Delayed Rejection . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Integration of Line Hypotheses . . . . . . . . . . . . . . . . 3.4.1 Generation of K Line Hypotheses . . . . . . . . . . . . . . . 3.4.2 Combination of Line Hypotheses into a Probability Map . . 3.4.3 Curvilinear Structure Extraction from Reduced Sampling Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . S AMPLING techniques with geometric priors have been exploited to detect multiple line segments in a scene [Lacoste 05, Tu 06, Schlecht 07, Batool 12, Predoehl 13, Jeong 14]. The marked point process (MPP) framework [

  Figure 3.1: A spatial point process is a mathematical model to analyze spatial point patterns using a finite number of parameters. More precisely, the spatial point process is a random variable which maps an unordered set of points in a finite space. We can compare different point patterns according to the randomness in the space.

  Figure 3.2: Comparison with the realization of point patterns obtain by (a) homogeneous Poisson process, (b) inhomogeneous Poisson process, and (c) repulsive Markov point process. For a finite space, homogeneous Poisson process employs a constant intensity measure, whereas inhomogeneous Poisson process employs intensity measure as a function. Both Poisson processes ignore correlation of points in the space. On the other hand, Markov point process constrains the local interaction of neighboring points using a non-negative density function. Here, repulsive constraint prevents overlapping between points.

  Figure 3.3: Line segments (s m , s n ) interact if they are close enough.

  10) where K is the total number of training data, b represents the bias of the hyperplane, and ξ ≥ 0 denotes slack variables for regression. The constant value C controls the relative importance of regularization ( 1 2 w w) and regression ( C K K i=1 ξ i ) terms (please refer to Chapter 4 for more detailed description of SVM).

  Figure 3.4: For the input image (a), we compare the manually labeled result (b), outcome of the binary SVM classifier (c) and the posterior probabilities (d). The proposed stochastic model (e) refines the weak points of the SVM classifier by sampling line segments with geometry prior constraints.

Figure 3 . 5 :

 35 Figure 3.5: For coupling feature c mn We measure (a) the distance of adjacent line segments d mn and (b) the angle differnece between line segments θ mn .

Figure 3 . 6 :

 36 Figure 3.6: To prevent congestion of line segments, we compute the overlapping area of pairwise line segments. We take into account either pixels on the line segment s m and its surroundings A(s m ), which are denoted by gray color in this figure. The criteria function Υ(s m , s m ) computes the number of pixels falling into the overlapping area of two line segments. We allow small overlapping between two line segments (a), whereas we reject line configuration (b).

  .7). We simulate a discrete Markov chain using a transition kernel K(s a , s b ) = Pr(X b = s b |X a = s a ) which updates state from X a to X b . To obtain an unique stationary probability p(s b |I), the Markov chain shall satisfy ergodicity [Chatelain 12]: p(s b |I) = sa∈Ψ p(s a |I)K(s a , s b ), ∀s b ∈ Ψ. (3.20) Thus, the Markov chain eventually reaches an equilibrium state which maximizes the proposed density in the configuration space. In practice, this property (3.20) can be implied by the detailed balance condition [Green 95, Chatelain 12] for any s a = s b ∈ Ψ:

  The sub-transition kernel is separated into a proposal probability Q m (s b |s a ) for moving from s a to s b and the corresponding acceptance ratio α m (s a , s b ). That is K m (s a , s b ) = α m (s a , s b )Q m (s b |s a ). With holding the detailed balance condition (3.21), we compute the acceptance ratio α m (s a , s b ) as follows [Hastings 70]:

  µ M (•) is Lebesgue measure for marked space M. That is µ M (M) = ( maxmin )(θ max -θ min ).

Figure 3

 3 Figure 3.8: The MCMC sampler perturbs the current configuration s via subtransition kernels: (a) birth kernel, (b) death kernel, and (c) linear-transform kernel. Schematic diagrams illustrate how the sub-transition kernels change topology of the corresponding configuration, and a state of the Markov chain is gradually updated in order to maximize the proposed probability density function. If a new configuration disagrees with the probability density function, the chain will remain at the previous state as shown in (d).

  Figure 3.9: Given configuration (a), if a line segment s proposed by LT kernel is rejected (b), the delayed rejection kernel searches for the nearest extremes in the rest of line segments (c). An alternative line segment s, which enforces connectivity, will be proposed by interpolation of the retrieved points (d).

Figure 3

 3 Figure 3.10: Given the input image (a), we compute the gradient magnitude (b). Mathematical morphology operator, path opening [Talbot 07], is applied on such gradient magnitude image (c). Linearity score of each pixel is drawn by the supervised feature learning algorithm [Becker 13] (d). We provide line hypotheses (e)-(g) associated with different hyperparameter vectors. Composition result (h) is equivalent to mixture probability density, and it highlights pixels corresponding to linear structures.

EE

  data (s i ) + sm∼sn E prior (s m , s n ; ω k ) 3: Compute a probability map P k via (3.29) 4: end for 5: Integrate K line hypotheses P S = 1 K K k=1 P k 6: Find the most promising hyperparameter ω k: k = argmax k={1,...,K} CC(P S , P k ) 7: Update the data likelihood energy (3.33) 8: Simulate Markov chain of the modified cost function over the reduced sample space S with ω k: data (s i ) + sm∼sn E prior (s m , s n ; ω k).

86 Figure 3 . 11 :

 86311 Figure 3.11: We provide intermediate sampling processes when the temperature parameter T t is decreasing. The results shown in first row are obtained without specifying seed segment. For the second row, we randomly set 20 seed segments and run the algorithm. For the third row, we initialize 20 line segments which are highly corresponding to underlying curvilinear structures. The algorithm converges toward almost the same solution regardless of the initial state.

Figure 3

 3 Figure 3.12:Precision-and-recall curves for pixelwise segmentation of curvilinear structures using path opening operator [Talbot 07] with different setups of length, supervised feature learning[Becker 13], baseline MPP, and the proposed method.

Figure 3

 3 Figure 3.13: We visualize the localization of the curvilinear structures on a Retina image. Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to ground truth (b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground truth.

Figure 3

 3 Figure 3.14: We visualize the localization of the curvilinear structures on a DNA image. Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to ground turth(b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground truth.

Figure 3 .

 3 Figure 3.15: We visualize the localization of the curvilinear structures on a wrinkle image. Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to ground truth (b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground-truth.

Figure 3

 3 Figure 3.16: We visualize the localization of the curvilinear structures on a Crack image. Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to ground truth(b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground truth.

Figure 4 . 1 :

 41 Figure 4.1: Comparing with the segmentation ((c), [Becker 13]) and the centerline detection ((d), [Sironi 14]) methods, the proposed algorithm is able to organize the topological importance of the curvilinear structure. The output rankings (e) highlight areas according to the geometric similarity of line segments and the underlying curvilinear structure. Our graphical model represents different levels of detail of the latent curvilinear structure with the minimum number of pixels. The proposed algorithm reconstructs the coarsest structure (f), and then iteratively adds fine branches (g). The final result (h) is obtained by a few iterations.

Figure 4 . 3 :

 43 Figure 4.3: Even a single outlier (see black dot) can dramatically change the characteristic of the hyperplane.

  Figure4.4: For the given linearly non-separable data points, we employ soft margin SVM which allows small classification errors around margin. A slack variable measures the such violation on the marginal areas. We modify the objective function of SVM to take into account the amount of margin violation with a parameter C. The hyperplane on the left figure is obtained using a small C so that it is less sensitive for classification errors and margin violations. On the other hand, the hyperplane on the right figure is more strict with respect to the margin violation by employing a larger C value.

Figure 4

 4 Figure 4.5: Overview of the curvilinear feature descriptor generation (Section 4.3) and ranking learning system (Section 4.4)

Figure 4 . 7 :

 47 Figure 4.7: Steerable filter banks f θ to compute curvilinear feature maps

Figure 4 . 9 :

 49 Figure 4.9: The topology of the spatial feature grouping is depicted in (a), where the corresponding weighting factors ω j , j ∈ N i of the neighboring line segments are denoted in blue colored numbers. The spatial grouping of the curvilinear features is able to improve the local consistency of the output rankings. Given image (b) is corrupted by irregular background textures. We visualize the output rankings: (c) is obtained without spatial feature grouping, whereas (d) shows the output rankings with the proposed spatial feature grouping.

0

  Figure 4.10: Example of the evaluation of the given line segment s i and corresponding ground truth patch g i . Blue and red pixels show superfluous and missed areas, respectively, whereas green pixels paint the agreed area between line segment and ground truth.

)

  Many options, such as correlation coefficient[Jeong 15a], Jaccard index [Jaccard 12], and Hamming distance[Hamming 50], are available to compute the function score.

  , z; w) = w Ψ(s, z). (4.19)

Figure 4 . 11 :

 411 Figure 4.11: Binary segmentation maps are built depending on statistics of the output rankings.

  (a) Binary segmentation map (b) Dissimilarity score map

Figure 4 .Algorithm 5

 45 Figure4.12: A binary segmentation map is created by thresholding output rankings at a deterministic value, where the value is obtained by the statistics of pixels on the curvilinear structure from the training dataset. While the structured output rankings contain the shape information of the line segments, we cannot exploit such useful information based on the binary segmentation map. The dissimilarity score map is generated by re-aggregation of the binary segmentation map. For each pixel on the binary segmentation map, we introduce the shape information which is already embedded in the output rankings. Then, we compute dissimilarity scores of the disassembled line segments and the binary segmentation map. Finally, we fill areas corresponding to the line segments with the dissimilarity scores.

Figure 4

 4 Figure 4.13: Toy example of the proposed curvilinear structure inference algorithm: (a) input image contains a curvilinear structure which is denoted by gray color; (b) subgraph G is induced from the structured output rankings; (c) and (d) show the intermediate processes of the 2-sweep algorithm starting from node t to find a diameter of the subgraph; (e) we assign 0 weight for all edges on the path; and (f) we repeat the process and add branches if the path length is larger than pre-defined length min .

Figure 4 .

 4 Figure 4.14: Intermediate steps of the curvilinear structure reconstruction for a retina image. We iteratively reconstruct the curvilinear structure according to topological importance orders. As the iteration goes on, detail structures (layers) appear.

  Figure 4.20. Moreover, Table 4.1 summarizes the Mean Average Precision (MAP) scores [Philbin 07, Jégou 08, Arandjelović 11] corresponding to Figure 4.20 for each dataset. We plot the curves by adjustment of the threshold value to represent the curvilinear structure. The dissimilarity score map, which is an intermediate result of the proposed algorithm, is used to plot Figure 4.20 and to compute Table

Figure 4

 4 Figure 4.15: Intermediate steps of the curvilinear structure reconstruction for a road network image.

Figure 4 .

 4 Figure 4.16: Experimental results on DRIVE dataset

Figure 4 .

 4 Figure 4.17: Experimental results on RecA dataset

Figure 4 .

 4 Figure 4.18: Experimental results on Aerial dataset

Figure 4 .

 4 Figure 4.19: Experimental results on Crack dataset

Figure 4 .

 4 Figure 4.20: Precision-and-recall curve of the curvilinear structure segmentation algorithms for each dataset. We plot the curve by controlling threshold values for each algorithm: [Frangi 98], [Law 08], [Becker 13], [Sironi 14] and dissimilarity scores of the proposed algorithm.

Figure 4 .

 4 Figure 4.21: We analyze the relationship between the average proportion of pixels to draw the curvilinear structure and the corresponding F 1 score for each dataset.

Figure 5 .

 5 Figure 5.22: Failure case of our stochastic model proposed in Chapter 3. Specifically, prior constraints penalize a large curvature of the local configuration, so that yield conflict at junction areas.

Figure A. 5 :

 5 Figure A.5: Example images of Crack dataset

Figure B. 2 :

 2 Figure B.2: Cas de défaillance du modèle stochastique proposé au chapitre 3. Plus précisément, les contraintes antérieures pénalisent une grande courbure de la configuration locale, pour reduire les conflicts dans les zones de jointure.

  Curvilinear structures are observed in different topological features under the various background textures. Left to right images show the latent curvilinear structure in blood vessels, DNA filaments, facial wrinkles, and road cracks. . . . . . . . . . . . . . . . . . . . . . . . . . . .

	List of Figures
	1.1
	h(•)	Hypothesis, decision function
	D	Training dataset for supervised learning
	φ(•)	Feature generating function
	Ψ(•, •)	Joint feature generating function
	z i	Latent variable
	ẑ	Estimated value obtained via hypothesis
	∆(z, ẑ)	Loss function to evaluate prediction
	ξ	Slack variable

2.1 The local orientation of the given image I is analyzed by performing image convolution with rotatable filter banks. Filtering responses accentuate pixels according to the specified orientation. The associated convolution filters F θ are depicted in left bottom of each image. . . . . 2.2 A candidate kernel is generated by the linear combination of shifted anisotropic Gaussian kernels. This approach is employed to analyze the size of underlying curvilinear structure intuitively. . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . Example of the feature maps of the morphological profile φ M with different length parameters L ∈ {50, 100, 150, 200} and those of the steerable filtering responses φ S with variable orientations θ ∈ {0• , 45 • , 90 • , 135 • }. The model parameters, which define the relative weights of those feature maps, are obtained via the ranking learning system. The higher ranking scores (denoted by red color) correspond to the latent curvilinear structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Steerable filter banks f θ to compute curvilinear feature maps . . . . . . 4.8 The baseline steerable feature vector φ S (s i ) is created by steerable filtering responses (F θ * I) as the increasing order of orientation. To take into account the direction of line segment s .13 Toy example of the proposed curvilinear structure inference algorithm: (a) input image contains a curvilinear structure which is denoted by gray color; (b) subgraph G is induced from the structured output rankings; (c) and (d) show the intermediate processes of the 2-sweep algorithm starting from node t to find a diameter of the subgraph; (e) we assign 0 weight for all edges on the path; and (f) we repeat the process and add branches if the path length is larger than pre-defined length min . . . . 4.14 Intermediate steps of the curvilinear structure reconstruction for a retina image. We iteratively reconstruct the curvilinear structure according to topological importance orders. As the iteration goes on, detail structures (layers) appear. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Intermediate steps of the curvilinear structure reconstruction for a road network image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Experimental results on DRIVE dataset . . . . . . . . . . . . . . . . . 4.17 Experimental results on RecA dataset . . . . . . . . . . . . . . . . . . 4.18 Experimental results on Aerial dataset . . . . . . . . . . . . . . . . . . 4.19 Experimental results on Crack dataset . . . . . . . . . . . . . . . . . . Cas de défaillance du modèle stochastique proposé au chapitre 3. Plus précisément, les contraintes antérieures pénalisent une grande courbure de la configuration locale, pour reduire les conflicts dans les zones de jointure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

4.5 Overview of the curvilinear feature descriptor generation (Section 4.3) and ranking learning system (Section 4.4) . . . . . . . . . . . . . . . . 4.6 i , we multiply permutation matrix P θ i by the baseline φ S (s i ). . . . . . . . . . . . . . . . . . . . . 4.9 The topology of the spatial feature grouping is depicted in (a), where the corresponding weighting factors ω j , j ∈ N i of the neighboring line segments are denoted in blue colored numbers. The spatial grouping of the curvilinear features is able to improve the local consistency of the output rankings. Given image (b) is corrupted by irregular background textures. We visualize the output rankings: (c) is obtained without spatial feature grouping, whereas (d) shows the output rankings with the proposed spatial feature grouping. . . . . . . . . . . . . . . . . . . 4.10 Example of the evaluation of the given line segment s i and corresponding ground truth patch g i . Blue and red pixels show superfluous and missed areas, respectively, whereas green pixels paint the agreed area between line segment and ground truth. . . . . . . . . . . . . . . . . . . . . . . 4.11 Binary segmentation maps are built depending on statistics of the output rankings. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 A binary segmentation map is created by thresholding output rankings at a deterministic value, where the value is obtained by the statistics of pixels on the curvilinear structure from the training dataset. While the structured output rankings contain the shape information of the line segments, we cannot exploit such useful information based on the binary segmentation map. The dissimilarity score map is generated by re-aggregation of the binary segmentation map. For each pixel on the binary segmentation map, we introduce the shape information which is already embedded in the output rankings. Then, we compute dissimilarity scores of the disassembled line segments and the binary segmentation map. Finally, we fill areas corresponding to the line segments with the dissimilarity scores. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44.20 Precision-and-recall curve of the curvilinear structure segmentation algorithms for each dataset. We plot the curve by controlling threshold values for each algorithm: [Frangi 98], [Law 08], [Becker 13], [Sironi 14] and dissimilarity scores of the proposed algorithm. . . . . . . . . . . . 4.21 We analyze the relationship between the average proportion of pixels to draw the curvilinear structure and the corresponding F 1 score for each dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.22 Failure case of our stochastic model proposed in Chapter 3. Specifically, prior constraints penalize a large curvature of the local configuration, so that yield conflict at junction areas. . . . . . . . . . . . . . . . . . . . A.1 Example images of DNA dataset . . . . . . . . . . . . . . . . . . . . . A.2 Example images of Drive dataset . . . . . . . . . . . . . . . . . . . . A.3 Example images of Wrinkle dataset . . . . . . . . . . . . . . . . . . A.4 Example images of Aerial dataset . . . . . . . . . . . . . . . . . . . A.5 Example images of Crack dataset . . . . . . . . . . . . . . . . . . . . B.1 Des structures curvilignes sont observées avec différentes caractéristiques topologiques et diverses textures de fond. De gauche à droite des images montrent une structure curviligne latente de vaisseaux sanguins, de filaments d'ADN, de rides du visage et de fissures sur une route. . . . B.2

  • S.-G. Jeong, Y. Tarabalka, N. Nisse and J. Zerubia, Inference of curvilinear structure based on learning a ranking function and graph theory, IJCV (under review).

Chapter2 Contents 2.1 Curvilinear Feature Descriptors . . . . . . . . . . . . . . . . 14

	• S.-G. Jeong, C. Lee and C.-S. Kim, Exemplar-based frame rate up-conversion
	with congruent segmentation, in Proc. IEEE ICIP, pp. 845-848, September
	2012.
	• S.-G. Jeong and C.-S. Kim, Feature-preserving thumbnail generation based on
	graph cuts, in Proc. IEEE ICIP, pp. 1081-1084, September 2011.
	2.1.1 Image gradient profiles . . . . . . . . . . . . . . . . . . . . . . 15
	2.1.2 Image morphological profiles . . . . . . . . . . . . . . . . . . 17
	2.2
	• S.-G. Jeong, Y. Tarabalka, N. Nisse and J. Zerubia, Progressive curvilinear
	structure reconstruction with structured ranking learning and graph algorithm,
	in Proc. IEEE CVPR, June 2016 (in preparation).

Existing Curvilinear Structure Reconstruction Methods . 20

	2.2.1 Curvilinear structure segmentation . . . . . . . . . . . . . . . 21
	2.2.2 Graph-based models . . . . . . . . . . . . . . . . . . . . . . . 24
	2.2.3 Stochastic image models . . . . . . . . . . . . . . . . . . . . . 25
	2.2.4 Machine learning based approches . . . . . . . . . . . . . . . 26

2.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

•

  Detection of the defects in the asphalts. Blood vessel extraction algorithms have been proposed to aid diagnosis of vascular diseases in retinal images [Chaudhuri 89, Frangi 98, Zana 01, Jiang 03, Staal 04, Al-Diri 09, Rigamonti 12]. Retinal images are scanned by ophthalmoscope, which is an optical microscopic device. As the brightness of the scanned images varies according to the illumination environment, image processing techniques are employed to enhance the visibility of blood vessels on the retinal images. The curvilinear structures of vascular images are mainly composed of diameter-varying vessels with irregular shapes. Curvilinear structure segmentation algorithms exploit low-level image features and multiscale manners to analyze vascular networks [Frangi 98, Law 08].

	Moreover, curvilinear objects, such as dendrites of neurons and RecA-DNA fila-
	ments, are observed in biological images. Neuroimage scientists have attempted to
	analyze such curvilinear structures in microscopic images [

Table 4 .

 4 1: Mean Average Precision (MAP) scores of each dataset for quantitative evaluation. Boldfaced numbers are used to indicate the best performance in each test.

		DRIVE RecA Aerial Cracks
	[Frangi 98]	0.61	0.55	0.38	0.074
	[Law 08]	0.73	0.33	0.25	0.14
	[Becker 13]	0.83	0.70	0.88	0.37
	[Sironi 14]	0.68	0.64	0.89	0.36
	Proposed	0.67	0.77	0.76	0.25

Table 4 .

 4

		DRIVE	RecA	Aerial	Cracks
	[Frangi 98] 0.33 / 37.09% 0.33 / 16.14% 0.32 / 22.19% 0.056 / 28.64%
	[Law 08]	0.43 / 19.60% 0.21 / 29.82% 0.25 / 29.51% 0.085 / 33.00%
	[Becker 13] 0.50 / 12.67% 0.45 / 9.39% 0.53 / 9.76%	0.23 / 32.13%
	[Sironi 14]	0.55 / 5.41% 0.50 / 5.47% 0.55 / 0.83%	0.27 / 17.12%
	[Jeong 15a] 0.59 / 0.97% 0.57 / 0.35% 0.47 / 0.46%	0.12 / 0.34%
	Proposed	0.36 / 1.07% 0.59 / 0.35% 0.59 / 0.57% 0.38 / 0.17%

2: F 1 scores and the average proportion of the pixels being a part of curvilinear structure (%) are computed. Boldfaced numbers are used to show the best F 1 score in each test. Gray colored cells in this table denote the algorithm using the minimum number of pixels to represent the curvilinear structure.

  based on machine learning. 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. Specifically, since RecA dataset has only four training images, there is a high risk of over-fitting to the training set. It is remarkable that the proposed algorithm shows good performance for all datasets without over-fitting problem when compared to other learning based algorithms[Becker 13, Sironi 14].

https://www.biochem.wisc.edu/faculty/inman/empics/dna-prot.htm

http://hugues.zahlt.info/91.html

http://cvlab.epfl.ch/page-108936-en.html

http://cvlab.epfl.ch/software/centerline-detection

Notre modèle stochastique proposé au chapitre 3 a souvent trouvé des segments déconnectés de jonctions en T (voir Figure B.2). Car le modéle est conçu pour rejeter de grandes courbures dans la configuration. Nous pouvons remédier à cet inconvénient en utilisant une autre couche pour représenter la structure curviligne sous une forme différente, par exemple, des réseaux de courbes de Bézier[Bama 15] ou une séquence de pixels. On pourrait définir une fonction, qui soit inversible, transformée à partir d'un groupe de segments de ligne en une représentation souhaitée. Nous supposons une forte connectivité des pixels appartenant au même groupe. Pour mettre en oeuvre ce modèle, nous devons concevoir un nouveau noyau de transition pour l'échantillonneur MCMC à saut réversible.En ce qui concerne notre second modèle, fondé sur l'apprentissage automatique, nous pourrions améliorer le descripteur de fonction pour tenir compte de la forme

Chapter3

Chapter 4

Inference of Curvilinear Structure based on Learning a Ranking Function and Graph Theory Abstract

To detect curvilinear structures in natural images, we propose a novel ranking learning system and an abstract curvilinear shape inference algorithm based on graph theory. We analyze the curvilinear structures as a set of small line segments. In this chapter, the rankings of the line segments are exploited to systematize the topological feature of the curvilinear structures. Structured Support Vector Machine is employed to learn the ranking function that predicts the correspondence of the given line segments and the latent curvilinear structures. We employ image gradient profiles and morphological profiles to describe curvilinear features. Also, we propose an orientation-aware feature descriptor and a feature grouping operator to improve the structural integrity during the learning process. To infer the curvilinear structure, we build a graph using the output rankings of the line segments. We progressively reconstruct the curvilinear structure by looking for paths between remote vertices in the graph. Experimental results show that the proposed algorithm faithfully detects the curvilinear structures within various datasets.

Conclusions and Perspectives

information regarding pairwise relation of line segments. Without this information, it was difficult to define a loss function which measures the quality of the estimated parameter for the cost function.

We can also consider to speed up the proposed algorithms by employing parallel programming. Specifically, the optimization process of our stochastic model supports the parallel implementation. The computational time can be further reduced by applying a parallel Monte Carlo sampler [Verdié 14].

Appendix

Appendix A Data Sets T HIS chapter describes the datasets used for the experiments conducted in this thesis. We evaluate the proposed curvilinear structure reconstruction models on various types of images, such as microscopic images of RecA-DNA filaments; retinal images of blood vessels; facial wrinkles; satellite images of road networks; and defects in the asphalt. The details regarding each dataset with the representative images are presented below.

A.1 DNA: Microscopic images of RecA-DNA filaments

In biological chemistry, negative stain of electron microscopy visualizes protein complexes [Jain 95]. We collect electron microscopic images of RecA proteins which contain filament structure (see

A.3 Wrinkle: Facial wrinkles

We collect facial wrinkle images of well known people on the Internet, such as actors and politicians. We also tested our methods on the facial skin images which are collected by world a leading dermatology company; however, these images cannot be presented in this thesis due to the permission issues. Facial wrinkles can be occluded by hairs and accessories or covered by makeup. Thus, we mainly test the algorithms to detect the wrinkles on the forehead part, where the corresponding areas are manually labeled before applying the algorithms.

A.4 Aerial: Satellite images of road networks

We test satellite images containing road networks to evaluate the curvilinear structure reconstruction algorithms. The dataset contains 14 remote sensing images of road networks and divided as 7 images for training and 7 images for the testing, respectively. Curvilinear structures associated with road network are partially occluded by the trees or cars on the road. This dataset is courtesy of the authors of [Sironi 14] and available at their website 2 . Some images of this dataset are provided in

Résumé étoffé de la thèse en Français B.1 Introduction

La recherche et les industries liées à la vision par ordinateur ont eu un développement rapide ces dernières années en raison de la fiabilité des capteurs et du développements de la puissance de calcul, cela en analysant différents types d'images. De l'imagerie microscopique à l'imagerie satellitaire, nous observons que des structures curvilignes existent sous différentes formes. Il est donc nécessaire de détecter les structures curvilignes latentes afin d'améliorer la visibilité des objets recherchés, d'extraire des informations biologiques, et d'analyser l'information géographique. Alors que l'homme peut intuitivement distinguer ces structures curvilignes latentes, la détection automatique de ces structures curvilignes est difficile à ce jour. Les difficultés sont dues au fait que ces structures curvilignes sont présentées sous différentes topologies. Ainsi, il est impossible de définir ces structures curvilignes complexes avec seulement quelques paramètres. Selon l'environnement de l'acquisition, les structures curvilignes peuvent être cachées sous les textures de fond avec un faible contraste. Par conséquent, la représentation des données et de nombreux algorithmes de vision par ordinateur, pour détecter les structures curvilignes, ont été développés pour des applications spécifiques. La figure B.1 montre des applications variées pour l'utilisation des structures curvilignes.

En vision bas niveau, des chercheurs [Freeman 91, Perona 95, Jacob 04] ont développé des filtres de convolution pour extraire les caractéristiques curvilignes de la texture du fond. L'idée principale derrière cette conception de filtre est de créer des modèles de forme de ligne pour extraire des informations de gradients orientés localement. Frangi et al. analysent les valeurs propres de la matrice de hessien pour une image l'image correspondante par des informations de gradient orienté et par filtrage morphologique d'image. Nous avons également réalisé que des a priori géométriques sont essentiels pour améliorer les performances en inférant des informations manquantes.

Nous avons proposé deux modèles de reconstruction de structures curvilignes différentes : un modèle stochastique dans le chapitre 3 et un modèle fondé sur l'apprentissage automatique dans le chapitre 4. Le chapitre 3 a commencé par proposer une densité de probabilité de segments de ligne d'images données. Afin de maximiser la densité proposée, nous avons utilisé l'échantillonneur RJMCMC avec un rejet retardé. Le modèle stochastique proposé nécessite le choix d'hyperparamètres pour décider de l'importance relative entre les contraintes géométriques. Au lieu d'estimer les valeurs exactes de ces hyperparamétres, nous avons simulé plusieurs chaînes de Markov pour obtenir des hypothèses de ligne. Un algorithme d'agrégation de ces hypothèses de ligne est proposé pour réduire l'espace d'échantillonnage et trouver l'ensemble des paramètres les plus prometteurs.

Le chapitre 4 a proposé une nouvelle méthode de représentation de la structure curviligne fondée sur la théorie des graphes qui exploite le classement de sortie des segments de ligne. L'algorithme d'inférence de structures curvilignes proposé visualise l'importance topologique en fonction de la longueur de la structure. Nous avons appris une fonction de classement afin de prédire la correspondance entre les segments de ligne et la structure curviligne latente. L'algorithme proposé reconstruit efficacement la structure curviligne avec un nombre minimum de pixels utilisés. Pour chaque chapitre, nous avons comparé les performances quantitatives et qualitatives de l'algorithme proposé avec des méthodes concurrentes de l'état de l'Art.

B.4 Perspectives