Alena Schmidt
email: alena.schmidt@ipi.uni-hannover.de

Florent Lafarge
email: florent.lafarge@inria.fr

Claus Brenner
email: claus.brenner@ikg.uni-hannover.de

Franz Rottensteiner
email: rottensteiner@ipi.uni-hannover.de

Christian Heipke
email: heipke@ipi.uni-hannover.de

Forest point processes for the automatic extraction of networks in raster data

Keywords: Spatial point processes, RJMCMC, graphs, digital terrain models

come

Forest point processes for the automatic extraction of networks in raster data

Introduction

The automatic extraction of networks from images or lidar data is of great interest in various disciplines. For instance, in remote sensing data, roads and rivers represent networks. In medical data, neurons and vessels can be described as network. Knowledge about the appearance of these objects helps to update road maps, to make predictions in case of flood events and to provide information for automated diagnostic systems.

For the automatic extraction of objects, different techniques from the field of image analysis are employed, frequently with the objective to integrate apriori knowledge into the model. The knowledge is often expressed in terms of probabilities and, thus, qualifies probabilistic strategies for object extraction.

In image analysis, a well-known approach is the use of Markov Random Fields and Conditional Random Fields as introduced by [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF] and [START_REF] Kumar | Discriminative Random Fields: A discriminative rramework for contextual interaction in classification[END_REF]. Here, the image is represented by a graph whose nodes correspond to the pixels or segments; the edges indicate local context. Knowledge about the objects in the image may be integrated by favoring similar classes for pixels in a local neighborhood [START_REF] Li | Markov random field modeling in computer vision[END_REF].

However, it is difficult to integrate more global constraints of the objects, e.g. concerning their shape.

In contrast, probabilistic model-based approaches express knowledge about the objects in a more holistic way. While this requires a careful modeling of the objects in order to stay general enough, it allows the integration of object characteristics beyond pixel-or segment-based relations. In this context, the method of spatial point processes has been shown to achieve good results in various disciplines, e.g. for object detection in remote sensing using raster data ( [START_REF] Benedek | Moving target analysis in ISAR image sequences with a multiframe marked point process model[END_REF], [START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF], [START_REF] Lafarge | Geometric feature extraction by a multimarked point process[END_REF], [START_REF] Tournaire | An efficient stochastic approach for building footprint extraction from digital elevation models[END_REF], [START_REF] Ortner | Building outline extraction from digital elevation models using marked point processes[END_REF]), point clouds ( [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF], [START_REF] Huang | A generative statistical approach to automatic 3D building roof reconstruction from laser scanning data[END_REF], [START_REF] Descombes | Marked point process in image analysis[END_REF]) or for the object detection in terrestrial images ( [START_REF] Ge | Marked point processes for crowd counting[END_REF], [START_REF] Ripperda | Reconstruction of façade structures using a formal grammar and RjMCMC[END_REF]).

Such spatial point processes benefit from (1) their flexibility in integrating knowledge about the objects and their relation to other objects, (2) their variability of the number of objects which is not restricted in the sampling process, and (3) the determination of a global optimal object configuration.

In this paper, we aim to take advantage of these benefits in order to detect networks in raster data. The paper is structured as follows: First, state-of-theart approaches dealing with network detection are reviewed (Section 1.1). Then, we describe the mathematical foundation of stochastic optimization with spatial point processes using a Reversible jump Markov chain Monte Carlo sampler coupled with simulated annealing (Section 2). In Section 3, we present the models employed for detecting networks of rivers and tidal channels in digital terrain models (DTMs), but beyond the detection of different types of networks in optical images. Our experimental setting and the results for different data sets are described in Sections 4 and 5, respectively. Finally, conclusions and perspectives for future work are presented in Section 6.

Related work

Several methods have been developed for the automatic extraction of networks in raster data. In the following, we focus on methods for the detection of (1) rivers and tidal channels -using 3D input data and integrating physical knowledge about the flow direction of water and (2) roads in remote sensing and vessels in medicine -using 2D images and integrating knowledge about the gray values and characteristics of the network.

River extraction. For some hydrological tasks, it is necessary to describe the flow of the water in the whole scene, e.g. in order to enable flood predictions. This is typically done with flow routing algorithms using DTMs in the form of raster data as input data, see [START_REF] Gruber | Land-surface parameters and objects in hydrology[END_REF] and [START_REF] Wilson | Water in the landscape: a review of contemporary flow routing algorithms[END_REF]. Generally, these algorithms calculate the flow of water in the scene which is related to the catchment area of each pixel in the DTM. Considering that the movement is mainly driven by gravity and neglecting other influences such as the property of the materials and surface roughness, the flow of the water is only allowed to neighboring pixels with a lower height within these approaches. Two groups of flow routing algorithms can be distinguished: (1) single flow methods (e.g. D8 [START_REF] O'callaghan | The extraction of drainage networks from digital elevation data[END_REF], R8 [START_REF] Fairfield | Drainage network from grid digital elevation models[END_REF], Kinematic Flow [START_REF] Lea | An aspect-driven kinematic routing algorithm[END_REF]) where the movement is restricted to one neighboring cell and (2) multiple flow algorithms (e.g. TOPMODEL [START_REF] Quinn | The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models[END_REF], FD8 [START_REF] Freeman | Calculating catchment area with divergent flow based on a regular grid[END_REF], MFD [START_REF] Holmgren | Multiple flow direction algorithms for runoff modelling in grid based elevation models: An empirical evaluation[END_REF]) where the flow of the water can be subdivided into multiple neighboring pixels with a lower height. An advantage of these methods is the computational efficiency; some of them are implemented in GIS software packages. However, they focus on steep terrain and are not fully transferable to hydrological applications with nearly horizontal terrain such as tidal channel networks as shown in [START_REF] Lohani | Application of airborne scanning laser altimetry to the study of tidal channel geomorphology[END_REF].

The characteristics of tidal channels -structures similar to rivers in Wadden Sea areas -differ from those of normal rivers, because the flow direction changes four times a day due to the tides. As a consequence, specific methods for the automatic detection of these networks have been developed, mainly based on image analysis. For instance, [START_REF] Fagherazzi | Tidal networks -Automatic network extraction and preliminary scaling features from digital terrain maps[END_REF] consider the height and the curvature of the terrain in DTM at different scales and combine the results in a thresholdbased approach. The method detects most of the channels in the input data, but fails for some of the small channels and does not generate a completely connected network. In contrast, [START_REF] Mason | Extraction of tidal channel networks from airborne scanning laser altimetry[END_REF] develop an approach starting with low level image processing operations such as edge detectors which are subsequently processed in order to find a channel network. For that purpose, parallel edges corresponding to both channel borders are searched and median axes of the channels are determined and linked based on their directions. The authors also combine the approach with optical input data which, however, does not improve the results [START_REF] Lohani | Extraction of tidal channel networks from aerial photographs alone and combined with laser altimetry[END_REF].

Network extraction with local strategies. A far larger number of network extraction approaches can be found in the field of road detection, see e.g. [START_REF] Mayer | A test of automatic road extraction approaches[END_REF], and the diagnosis of networks and trees in medical data such as neurons or blood vessels, see [START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF] and [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF]. Here, optical images or 3D image stacks are used as input data. Some approaches solely consider local characteristics of the network such as geometric and radiometric features in a local neighborhood. These methods benefit from computational efficiency, but have the disadvantage that they are not robust against noise and fail in case of occlusions of the objects in the images. Among this group of approaches, some methods can be characterized as tracking approaches starting from some seed points and directions and expanding the network by iteratively adding points and paths. This is done using a Kalman filter strategy by [START_REF] Vosselman | Road tracing by profile matching and Kalman filtering[END_REF] and [START_REF] Movaghati | Road extraction from satellite images using particle filtering and extended Kalman filtering[END_REF] or by defining statistical tests concerning the width, direction and curvature of roads [START_REF] Geman | An active testing model for tracking roads in satellite images[END_REF]. In [START_REF] Hu | Road network extraction and intersection detection from aerial images by tracking road footprints[END_REF], segments are gradually added to the network by analyzing their shape and deriving the dominant directions of streets. For medical data tracking [START_REF] Chothani | Automated tracing of neurites from light microscopy stacks of images[END_REF] and [START_REF] Bas | Principal Curves as Skeletons of Tubular Objects[END_REF] analyze features based on the Hessian matrix of the gray values in a local neighborhood or the straightness of detected lines in the network. If many streets or vessels are close to each other, a disadvantage of tracking approaches is that the network may be expanded in an incorrect way.

Many approaches define a set of rules for delineating road or vessel networks.

Their general strategy is to detect parts of the network starting from an image segmentation ( [START_REF] Grote | Road network extraction in suburban areas[END_REF], [START_REF] Gamba | Improving urban road extraction in high-resolution images exploiting directional filtering, perceptual grouping, and simple topological concepts[END_REF], [START_REF] Poullis | Delineation and geometric modeling of road networks[END_REF], [START_REF] Mena | An automatic method for road extraction in rural and semi-urban areas starting from high resolution satellite imagery[END_REF]) or image classification ( [START_REF] Mnih | Learning to detect roads in high-resolution aerial images[END_REF], [START_REF] Marín | A new supervised method for blood vessel segmentation in retinal images by using graylevel and moment invariants-based features[END_REF], [START_REF] Zhang | Towards an operational system for automated updating of road databases by integration of imagery and geodata[END_REF]) and to group these parts based on prior knowledge. In general, the knowledge is integrated in a heuristic way and requires the tuning of a large set of parameters for each scene.

Network extraction with global strategies. In contrast to local methods, global approaches evaluate paths between seed points in the entire input data set and optimize the network for the whole scene. Active parametric contour models (snakes) [START_REF] Kass | Snakes: Active contour models[END_REF] first initialize a contour representing the network. This contour is deformed by internal forces (describing knowledge about the smoothness of the contour) and external forces (constraining the network to the data). For road networks the internal and external forces are described by constraints about the linearity and width of streets or about radiometric features in optical images [START_REF] Gruen | Road extraction from aerial and satellite images by dynamic programming[END_REF] or SAR data [START_REF] Bentabet | Road vectors update using SAR imagery: A snake-based method[END_REF]. [START_REF] Wang | Novel 4-D open-curve active contour and curve completion approach for automated tree structure extraction[END_REF] also allow the initialized network to expand by adding pixels in a local neighborhood in the case of similar features of these pixels. There are different possibilities to extend the parametrization of the contours, for instance by considering the width of the roads with Ribbon snakes [START_REF] Baumgartner | Automatic road extraction based on multi-scale, grouping, and context[END_REF] or modeling the contour by a graph in network snakes [START_REF] Butenuth | Network snakes: Graph-based object delineation with active contour models[END_REF] and, thus, benefit from the direct consideration of the network topology. However, the modeling by snakes requires a good initialization of the network and, thus, makes these methods more suitable to applications where an approximate contour is given or extracted in a preprocessing step.

Another group of global methods is based on minimum cost paths. By using graph structures, the topology is modeled explicitly. For road networks, one of the first approaches using minimum cost is the method of [START_REF] Fischler | Detection of roads and linear structure in low resolution aerial imagery using a multi-source knowledge integration technique[END_REF]. The results of different line detectors are combined and a score function indicating the probability belonging to the network is evaluated. The A* algorithm is employed to find the optimal path, i.e. the one having the highest scores (equivalent to the minimum costs). The minimum cost path is also used by [START_REF] Gerke | Graph-supported verification of road databases[END_REF] and [START_REF] Hinz | Automatic extraction of urban road networks from multi-view aerial imagery[END_REF] in order to connect road segments from a pre-processing step or to find the best path in classified image data [START_REF] Türetken | Automated reconstruction of tree structures using path classifiers and mixed integer programming[END_REF]. By evaluating all pixels along the edges, the methods based on minimum cost path can successfully bypass single noise pixels in the data. Moreover, these methods provide not only the geometry, but also the topology of the network. However, in general, all seed pixels have to be added to the graph, which may result in topological errors. Another strategy of network detection is to embed the problem in an image classification using graphical models based on Markov Random Fields or Conditional Random Fields. [START_REF] Tupin | Detection of linear features in SAR images: Application to road network extraction[END_REF] and [START_REF] Katartzis | A model-based approach to the automatic extraction of linear features from airborne images[END_REF] first generate line segments by using different line detectors and topological filters. [START_REF] Montoya-Zegarra | Mind the gap: modeling local and global context in (road) networks[END_REF] or [START_REF] Wegner | Road networks as collections of minimum cost paths[END_REF] calculate probabilities for each pixel or superpixel, respectively, for being part of the road network and, then, randomly choose pairs of pixels with high probabilities and connect them by optimal paths. By using higher order clique potentials global context is integrated. However, global constraints of the objects, e.g. concerning their shape, are difficult to integrate in such an approach.

Network extraction by stochastic approaches. The network detection by spatial point processes is used for road network detection for instance by [START_REF] Stoica | A Gibbs point process for road extraction from remotely sensed images[END_REF]. The authors introduce a model consisting of line segments having a certain width.

Two line segments are considered to be connected if the smallest Euclidean distance between their endpoints is smaller than a predefined threshold. In the optimization process, the line segments are modified by changing their lengths and widths. Moreover, line segments are added to or removed from the configuration. For the evaluation of each object configuration, segments which are not connected to other segments on both sides are penalized. They are also penalized if they overlap, are too short, or enclose a small angle. The conformity of the object hypotheses with the input data is evaluated based on the homogeneity of the pixel gray values inside the segments and the gray value differences to the pixels outside the segment. [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF] refine the method and avoid constant penalization terms by evaluating the segment configuration depending on the angle and the distance between two adjacent segments. An expansion of the model of [START_REF] Stoica | A Gibbs point process for road extraction from remotely sensed images[END_REF] to 3D data is realized by [START_REF] Stoica | A three-dimensional object point process for detection of cosmic filaments[END_REF] and [START_REF] Kreher | Detektion von Hirnnervenfasern auf der Basis von diffusionsgewichteten Magnetresonanzdaten[END_REF], modeling the objects by cylinders instead of line segments.

The methods based on marked point processes mentioned above have in common that they require the definition of geometric constraints in order to regard segments as connected, and, thus, enforce the connectivity of the network. Moreover, they do not directly provide the topology of the network.

In contrast, this is done in [START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF] by modeling the network as a planar graph.

In the sampling process, the nodes representing junction-points are iteratively connected by edges, which correspond to line segments and may have different widths. The edges are evaluated based on the homogeneity of the gray values inside the segments and the gradient magnitudes on the segment borders. The authors penalize graph configurations with non-connected components, atypical intersection angles and atypical numbers of outgoing edges (both compared to training data). However, the applications are restricted to 2D image data, and the method does not integrate physical knowledge. Moreover, the graph in the approach of [START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF] contains cycles which does not seem appropriate as object model for our main target application: the extraction of river networks. In general, river networks are structured as tree, where two or more streams converge and form a bigger stream in the network.

Contributions

In this paper, we propose a new type of spatial point process for the detection of line networks. Similarly to [START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF], our spatial point process is designed to sample planar graphs. In our case, graphs are guaranteed to be acyclic so that the sampled structure is a tree or a set of trees. We call this process a forest point process. This specificity requires us to define adequate perturbation operators. Moreover, our model differs from the approach of [START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF] by the energy function which is adapted for the detection of river networks in our model.

Here, we integrate physical knowledge about the network from 3D input data.

In comparison to our previous work ( [START_REF] Schmidt | Network detection in raster data using marked point processes[END_REF], [START_REF] Schmidt | A graph based model for the detection of tidal channels using marked point process[END_REF]), we present an enhanced model which is characterized by a reduced number of parameters, an extension of the energy model and an additional kernel. Our contribution is that (1) we develop a new spatial point process for capturing tree structures in images, (2) we take advantage of the 3D input data and, thus, facilitate the integration of physical knowledge about river networks, (3) we develop a hierarchical optimization strategy and (4) we show the transferability to river networks in different types of terrain, but also to different networks in image data.

Theoretical background

Spatial point processes

Point processes belong to the group of stochastic processes, see [START_REF] Møller | Statistical inference and simulation for spatial point processes[END_REF] and [START_REF] Daley | An introduction to the theory of point processes: Volume I: Elementary theory and methods, 2nd Edition[END_REF]. A point process is a sequence of random variables X t that takes values in a state space X . In our application t represent the iteration in the sampling process. In object extraction, point processes are used to find the most probable configuration of objects in a scene given the data. A point is described by its location v i = (x i , y i ). In a marked point process, a set of marks m i , i.e. a multidimensional random variable describing an object of a certain type at position v i , is added to each point. If we characterize an object l i = (v i , m i ) by its location and mark, a marked point process can be thought of as a stochastic model of configurations of an unknown number n of such objects in a bounded region S (here: the DTM, thus the objects l i exist in X = R 2 ).

There are different types of point processes, each of them associated with a measure which assigns a numerical value to each suitable subset B of S and, thus, makes different realizations of point processes quantitatively comparable.

For the Poisson point process -one of the fundamental point processes -the measure is given by the expected number of objects in each subset. The Poisson point process assumes a complete randomness of the objects and models the probability p(n) for the number of objects by a discrete Poisson distribution

p(n) = λ(B) n e -λ(B) n! . (1)
In ( 1) the parameter λ -also called intensity parameter -corresponds to the expected value for the number of objects in the subset B. In practice, the assumption of complete randomness of the object distribution is often not justified, and more complex models are postulated instead. These models are described by a probability density function g with respect to a reference point process, which is usually defined as the Poisson point process.

In our approach, g is expressed through a Gibbs energy U under the form g ∝ exp -U . In this way, we integrate interactions between adjacent objects by using Gibbs point processes, which are also applied e.g. in ( [START_REF] Stoica | A Gibbs point process for road extraction from remotely sensed images[END_REF], [START_REF] Mallet | A marked point process for modeling lidar waveforms[END_REF]). The Gibbs energy can be modeled by the sum of a data energy U d (X) and a prior energy U p (X):

U (X) = βU d (X) + (1 -β)U p (X) (2)
where X ∈ Ω is the object configuration from the set of all configurations Ω and β ∈ [0, 1] describes the relative influence of both energy terms. The data energy U d (X) measures the consistency of the object configuration with the input data.

The energy U p (X) introduces prior knowledge about the object layout; our models for these two energy terms are described in Section 3.3. The optimal configuration X of objects can be determined by minimizing the Gibbs energy U (X), i.e. X = arg min X∈Ω U (X). This is typically done by coupling a Reversible jump Markov chain Monte Carlo (RJMCMC) sampler and a simulated annealing relaxation.

Reversible jump Markov chain Monte Carlo Sampling

225

Markov chain Monte Carlo (MCMC) methods belong to the group of sampling approaches. The special feature of MCMC methods is that the samples are not drawn independently, but each sample X t is drawn on the basis of a probability distribution that depends on the previous sample X t-1 . Thus, the sequence of samples forms a Markov chain in the space of possible configurations. If the number n of objects constituting the optimal configuration X were known and constant, it could be determined by MCMC sampling ( [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]).

RJMCMC is an extension of MCMC that can deal with an unknown number of objects and changes of the dimension of the parameters between two sampling steps [START_REF] Green | Reversible jump Markov Chain Monte Carlo computation and Bayesian model determination[END_REF]. The sampler proposes changes to the current configuration from a set of pre-defined types of changes. Each of the change types is associated with a density function Q i called proposition kernel. Each kernel Q i must be reversible, i.e. the inverse change must be possible [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF]. At each iteration t, a kernel Q i is chosen randomly according to a proposition probability p Qi which may depend on the kernel type. The configuration X is changed according to the kernel Q i , which results in a new configuration X . Subsequently, the Green ratio R ( [START_REF] Green | Reversible jump Markov Chain Monte Carlo computation and Bayesian model determination[END_REF])

is calculated: R = π(X ) π(X) p Q r i p Qi Q(X, h|X , h ) Q(X , h |X, h) det ∂f X→X (X, h) ∂(X, h) . (3)
In ( 4), p Qi and p Q r i are the probabilities for choosing kernel Q i and its inverse kernel Q r i , respectively. The kernel ratio Q(X,h|X ,h ) Q(X ,h |X,h) expresses the ratio of the probabilities for the change of the configuration from X to X and from X to X . In this context, possible changes in the parameter dimension are taken into account by introducing auxiliary dimension-matching variables h and h [START_REF] Green | Reversible jump Markov Chain Monte Carlo computation and Bayesian model determination[END_REF].

They depend on a transition function f which maps one configuration to another one, (X , h ) = f X→X (X, h). Its Jacobian determinant det ∂f X→X (X,h) ∂(X,h) is incorporated in equation 4. Both, f and h, have to be chosen depending on the application. π(X) is the target function that has to be optimized. By using the Gibbs energy in our application, equation 3 becomes

R = exp - U (X ) -U (X) T t p Q r i p Qi Q(X, h|X , h ) Q(X , h |X, h) det ∂f X→X (X, h) ∂(X, h) . (4)
We also coupled the RJMCMC sampler with simulated annealing. For that reason, the parameter T t referred to as temperature ( [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]) is introduced in equation ( 4). The sequence of temperatures T t tends to zero as t → ∞. Theoretically, convergence to the global optimum is guaranteed for all initial configurations X 0 if T t is reduced (cooled off) using a logarithmic scheme. In practice, a geometrical cooling scheme is generally introduced instead because it is faster and usually still gives a good approximate solution ([70], [START_REF] Tournaire | An efficient stochastic approach for building footprint extraction from digital elevation models[END_REF]).

Following [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] and [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] the new configuration X is accepted with an acceptance rate α and rejected with the probability 1 -α, where α is computed from

R using α = min(1, R). (5)
In practice, we sample a random number γ ∈ [0, 1]. If γ < α, the new configuration is accepted and X becomes the new configuration (X t+1 = X ). Otherwise the old configuration is maintained (X t+1 = X). The four steps of (1) choosing a proposition kernel Q i , (2) generating the new configuration X , (3) computing the acceptance rate α, and (4) accepting or rejecting the new configuration are repeated until a convergence criterion is achieved.

Extraction of line-networks

We use spatial point processes with RJMCMC sampling and simulated annealing as described in Section 2 to detect networks in raster data. First, we define the object representation of the networks (Section 3.1). Second, the types of change in the graph configuration which are applied during the sampling process are introduced (Section 3.2). Third, the energy function to be minimized during the global optimization is defined (Section 3.3).

Forest point process

Contrary to traditional marked point processes in which interactions between points are formulated in the energy through a neighboring relationship, we specify the adjacency of points in the space of marks directly. Each point of a configuration is thus associated to a set of point indices that gives the list of its neighboring points. This allows us to associate a unique adjacency graph to each configuration of points where the nodes of the graph correspond to the points and the edges correspond to pairs of neighboring points, as illustrated in Figure 1. As our target application is the extraction of rivers, we restrict the admissible topology to trees, i.e., undirected, acyclic graphs. Indeed we guarantee that water flows only in one direction by avoiding cycles in the graph.

Simultaneously, we accept that we cannot detect two channels circling round an island. Note that this particular case rarely occurs in our application. Since our algorithm is applied to scenes which may contain multiple independent rivers, the output is a collection of one or more disjoint trees, which is commonly known as a forest [START_REF] Bondy | Graph Theory[END_REF].

We define by forest point process a point process where each point is associated to a set of point indices so that the adjacency graph of any configuration of points is a forest. In other words, we restrict the form of the adjacency graph to undirected acyclic planar graphs to guarantee that any configuration of points can be represented as a set of trees.

In our approach, we enrich the space of marks so that each point is also associated to a parameter specifying the width of corresponding edges. A point l is thus defined by the set of variables l = (v, i 1 , ..., i k , w 1 , ..., w k ) where v = (x, y) refers to the position of the point, i 1 , ..., i k are the indices of its k adjacent points, and w 1 , ..., w k are the widths of the edges between the point and its adjacent points. When k = 1, the point is an extremity node. When k > 1, the point is a branch node that can be either associated to a change of direction (k = 2) or to a confluence of branches (k > 2), as illustrated in Figure 2.

Note that, when a pair of points are adjacent, they have an identical edge width in their set of marks. We denote by e the edge forming by two adjacent points in the adjacency graph.

Changes of the object configuration

In the sampling process, we modify the forest describing the network of objects in the raster data. For that purpose, potential changes (also referred to as perturbations) that can be applied as well as the corresponding kernels Q i must be defined. We allow three types of changes: (1) birth-and-death:

points are added to or removed from the current forest, (2) modification: the parameters of the points are modified, (3) split-and-merge: two points are split or merged which is done by merging the end nodes of two edges or by replacing their common end node by a new one for one of these edges while keeping the end node for the other edge unchanged. As mentioned in Chapter 2, each perturbation is required to be reversible, i.e. for each perturbation there has to exist an inverse perturbation; for instance, the death event is the inverse perturbation of the birth event. The corresponding kernel Q i indicates the probability for a specific change to a new configuration.

(

) 3
(2)

(1)

Figure 2: Types of nodes. The nodes correspond to either (1) extremities in the tree,

(2) changes of direction, or (3) confluences of branches.

Birth-and-death

In the optimization, nodes and edges are added to or removed from the current forest, which is accomplished by the birth-and-death kernel Q BD .

In the case of a birth event, the position of the new node (x j , y j ) is randomly chosen based on a probability map which assigns the probability for being part of the river network to each pixel in the DTM (see Section 3.4). Then, the node is added to the graph. In our model, each node has to be connected with at least one node in the forest. That is why we search for neighboring nodes in a local neighborhood r. We check if the connection to them does not result in a cycle in order to guarantee the acyclicity of the adjacency graph. Moreover, we check if the edges do not cross in the object space. From all neighboring nodes which fulfill these criteria, we randomly choose one and connect it to the new node by an edge (see Figure 3 top). If no neighboring points exist, the position of a second node (x (j+1) , y (j+1) ) is randomly chosen, again, based on the probability map and the second node as well as an edge connecting both nodes and, thus, a new tree are added to the forest (see Figure 3 bottom). Then, we randomly choose a width for the new edge. Here, we assume a uniform distribution for all values in the range of the minimal and maximal width w min and w max which are given for each scene. In total, one new edge and one or two new nodes are added to the forest in the birth event. An example is illustrated in Figure 3. The kernel ratio for this setting can be reduced to the ratio of the probabilities p(n)

for the number of nodes in the forest which is given by the Poisson distribution in equation 1. In this way, the Gibbs point process is set in reference to the Poisson point process. Then, the kernel ratio of the birth event is

Q BD (X, h|X , h ) Q BD (X , h |X, h) = p(n|n ) p(n |n) = λ n e -λ n ! • n! λ n e -λ =      λ n , n = n + 1 λ 2 n •(n -1) , n = n + 2. ( 6
)
In equation 6 parameters denoted with an apostrophe belong to the new configuration, parameters without an apostrophe to the old configuration, i.e. n is the number of nodes in the old configuration X and n is the number of nodes in the new configuration X . λ is the expected value for the number of nodes (see equation 1). The absolute value of the Jacobian determinant is 1 (see Appendix) and, thus, the Green ratio (equation 4) for the birth event for the case

n = n + 1 is R = exp - U (X ) -U (X) T t p QD p QB λ n (7)
where p QD and p QB are the probabilities for choosing the birth and the death 295 kernel, respectively.

In the case of a death event, we randomly choose a node with only one incident edge and remove it from the forest. We also delete its connection to the forest and remove the incident edge. If the adjacent node is not connected to a further node, we remove it, too. In total, one edge and one or two nodes 300 are removed from the forest in the death event.

Analogously to the birth event, only the probability for the number of nodes is required for the kernel ratio which is decremented by one or two in comparison to the old configuration. Thus, the kernel ratio for the death event is

Q BD (X, h|X , h ) Q BD (X , h |X, h) = λ n e -λ n ! • n! λ n e -λ =      n λ , n = n -1 n •(n -1) λ 2 , n = n -2. ( 8
)
Again, the Jacobian determinant is 1. Consequently, the Green ratio for the death event is given by

R = exp - U (X ) -U (X) T t p QB p QD n λ (9)
for n = n -1.

v 1 v 2 v 3 v 4 v 5 e 1 e 3 e 2 birth death v 1 v 2 v 3 v 4 v 5 v 6 e 1

Modification

For the modification of the parameters of the forest, we distinguish between three types of perturbation: We assume that all modifications are equally probable and randomly choose one of them. All types of modification are illustrated in Figure 4. For (1), we randomly generate a displacement vector within a local neighborhood. If the shift to the new position does not result in crossing edges, we change the node coordinates. For (2), we randomly change the width of an edge within a small range around the old width. For (3), we randomly choose a node and change its number of connections. This is done by adding or removing an edge to the forest. Both changes are equally probable, we randomly choose one of them. If we add a connectivity, similar to the birth event, we search for nodes in a local neighborhood r which can be connected to the chosen node without breaching the acyclicity or planarity of the adjacency graph. From all nodes which fulfill these criteria, we randomly choose one of them and add an edge connecting both nodes to the adjacency graph. If there does not exist a neighboring node, we do not change the configuration. For the reduction of the number of connections, we randomly choose a node and one of its incident edges. The edge is removed from the forest, if all nodes are still connected by another edge to the forest.

For the modification of the parameters, the Green ratio1 is

R = exp - U (X ) -U (X) T t p r QM o p QM o (10)
where p r QM o and p QM o are the probabilities for choosing the modification kernel and its reversible kernel, respectively. Note, that all coordinates and width changes in these events are set to be equally probable and, thus, the Kernel

ratio Q(X|X ) Q(X |X) is 1.

Split-and-merge

The third type of perturbation is a split-and-merge event where two nodes are split or merged. In the case of a merge event, a node is randomly chosen from the forest.

Then, nodes are searched in a local neighborhood. We check that merging two nodes does not lead to cycles or crossing edges in the forest. If more than one node fulfill these criteria, we randomly choose one of them. We set this node as new node for the edge and remove the chosen node (Fig. 5). The determinant of the Jacobian matrix for this setting is 1 and the kernel ratio is analogous to the birth event given by the ratio of the Poisson measure (equation 6). Then, the Green ratio of the split event is

R = exp - U (X ) -U (X) T t p QS p QM n λ . (11)
Conversely, in the case of a split event, a node with at least two incident edges is chosen randomly. Then, a local transition vector is randomly generated and a node of one of the incident edges is moved to the new position. A new node is added to the forest at this position. The Green ratio for this setting is

R = exp - U (X ) -U (X) T t p QM p QS λ n (12)
where p QM and p QS are the probabilities for choosing the merge kernel and the split kernel, respectively.

Energy model

For the evaluation of each configuration we model a Gibbs energy U (equation 2) which is composed of a data energy U d checking the consistency of the object configuration with the input data and a prior energy U p where we introduce prior knowledge about the configuration.

Data energy

The model of the data energy is motivated by our main target application:

the extraction of rivers and tidal channels in a DTM. These networks are characterized by locally lower heights than their surroundings. A typical profile of a tidal channel is given in Figure 6. Consequently, rivers and tidal channels are assumed to have high DTM gradient magnitudes on each bank and homogeneous heights in between. We model the data energy by

U d (X t ) = U g (X t ) + U h (X t ) (13)
where each term represents one of the two assumptions about rivers.

340

High gradient magnitudes at the borders: We determine the data energy from the DTM gradients at the edge borders by

U g (X t ) = ej ∈Xt   c 1 - 2 mj =1 1 n mj nj kj =1 ∇ ⊥ DT M k   . (14)
In [START_REF] Gruber | Land-surface parameters and objects in hydrology[END_REF],

∇ ⊥ DT M k = ∇ DT M k , m ⊥ j,k
is the component of the DTM gradient at boundary pixel k j in direction of the normal vector of the lateral boundary m j ∈ 1, 2 of the edge e j (see Fig. 7). The sum of the gradients is taken over the n j pixels k j along that border; all gradients have equal weights. We only take into account the two lateral borders m 1 , m 2 of the segment, potentially 345 corresponding to the river or channel banks. To ensure that the energy is only minimized, if the data term is above a predefined value, we introduce a constant c 1 .

Homogeneity of the heights: Perpendicular to the direction of the edge, we analyze the DTM heights by

m j = 3 m j = 1 m j = 2 m j = 4 k e j ∇ ⊥ DT M k ∇ DT M k
U h (X t ) = p h • ej ∈Xt max   0, -c 2 + 4 mj =3 σ z (m j )   , (15)
where σ z (m j ) is the standard deviations of the heights along the perpendicular borders m j ∈ 3, 4 of edge e j in the forest (see Fig. 7). p h is a weighting factor.

We introduce a constant c 2 in order to be less susceptible to noise in the data.

Since the data terms U g and U d are inconsistent with one another at the ends of the borders m j ∈ 3, 4 where U g forces the pixel to the pixels with the highest gradients and U d requires homogeneous heights of the pixels, we do not take into account 5 % of the pixels on both ends of the borders m j ∈ 3, 4 and calculate U h only for the 90 % of the pixels in the middle of the borders.

Note, that despite developed for river networks, the data energy is principally transferable to any network in raster data that consists of homogeneous gray values and has strong gradients at the borders. If the gray values of the network are higher than the gray values of the background, we use the complementary image as input data.

Prior energy

Prior knowledge is integrated into the model in order to favor certain object configurations independently of the input data. We model the prior energy by three terms

U p (X t ) = U o (X t ) + U s (X t ) + U f (X t ) (16)
where each term represents one of the characteristics of the network we intend to take into account in our model.

Non-overlapping segments:

We penalize an object configuration in which the edges overlap. In this way, the accumulation of objects in regions with high data energy can be avoided. For that purpose, we calculate the overlap area a of all combinations of edges e i and e j and their relative overlap area a a(ei) and a a(ej ) . The corresponding energy term is calculated by

U o (X t ) = p o i =j max a a(e i ) , a a(e j ) , (17)
whereby we only consider the larger relative overlap area of two edges e i and e j 365 and p o is a penalizing factor.

Graph connectivity: We want to obtain a configuration with only one tree and, thus, rate the connectivity of the nodes in the forest. For that purpose we penalize graph configuration consisting of two or more trees by

U s (X t ) = p s • (n sg -1) ( 18
)
where n sg is the number of trees and p s is a penalizing factor.

Physical consistency of the network: For river networks, we take into account the flow direction of water by analyzing the heights of the nodes and edges of the forest. The general idea is to take into account that water does not flow uphill. We integrate this physical knowledge using two criteria. First, each node needs to be connected to exactly one node with a lower height value but can have connections to an arbitrary number of nodes with a higher height value -this means that we do not consider deltas or islands in our model. Second, each pixel on the medial axis of an edge has to show the same trend in the data, i.e. from a node to its adjacent node all height values have to increase or decrease; otherwise the edge is penalized. Note that we allow small deviations σ from this trend in order to be less susceptible to noise in the data. Thus, an edge is not penalized if the height of each pixel follows the trend of the edge within a tolerance of σ. Combining both criteria, the prior energy is modeled by

U f (X t ) = ej ∈Xt p f • n f1 + n f2 (e j , σ) l(e j ) . (19)
In [START_REF] Quinn | The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models[END_REF], n f1 is the number of nodes in the forest violating the first criterion.

For each edge e j the number of pixels n f2 (e j ) violating the second criterion is counted and normalized by the length l(e j ) of the edge.

Note that we do not use an explicit shape prior in our energy. However, as illustrated in Fig. 1, tree branches are modeled by rectangles whose width and length live in a bounded domain. This constitutes an indirect way to impose some prior shape knowledge without adding an additional term in the energy which would require more delicate parameter tunings: When a perturbation operator modifies, adds or removes a tree branch, the perturbation is not fully random, but driven by some restriction w.r.t. the shape.

Hierarchical approach and probability map

We develop a strategy for a hierarchical extraction of networks. For that purpose, we apply the method on the input data of the raster data with different resolution steps. We then search for the edges with the largest width in the lowest resolution, keep the results, reduce the probability map for the detected pixels and complete the network in the higher resolution by detecting edges with lower widths. For river networks this means, that we first detect the biggest streams in the network and then search for the smaller streams.

The probability map indicates the probability to be a node to each pixel in the raster data. We calculate the probability map based on the threshold for simple raster data features, e.g. the heights in the DTM or the gray values in the image. For river networks, we also convolve the DTM with derivatives of the Gaussian and calculate the terrain curvatures based on the second derivative in the direction perpendicular to the gradients in the DTM. We apply a simple threshold method and assign a probability of 1 to all pixels below / above this threshold while the probability is set to 0.01 for the other pixels.

Experimental setting

Experimental setup

We tested our method on different data sets and for different applications.

First, the river and tidal channel network detection was evaluated. For that purpose, we applied our method for three DTMs. Synthetic data simulating an ideal tidal channel network were used in order to analyze the importance of each term in the energy function (Section 5.1). We also evaluated the influence of the generation of random numbers in the sampling process with this data sets (Section 5.2). The applicability to real data sets is shown based on two DTMs: a test site from the Wadden Sea with nearly horizontal terrain and a mountainous test site in Austria. On the one hand, we used this data for analyzing the influence of different parameter settings and, especially, for evaluating the influence of the integration of physical knowledge in the energy function (Section 5.3). On the other hand, we compared our results to those from a standard flow routing algorithm (Section 5.4). Finally, the transferability to different types of networks in images is shown (Section 5.5).

For the data sets of the Wadden Sea and Austria, we applied the hierarchical approach and searched for rivers in two resolution steps of the DTM. Moreover, for all experiments we used a geometrical cooling scheme in the simulating annealing process by reducing the temperature T t = T t0 • (d f ) t with a decreasing factor d f < 1. We also implemented a logarithmic cooling scheme, which theoretically guarantees to achieve the global optimum and test it for the synthetic data set. Note that we do not implement a stop criteria for the optimization, but rather set the number of iterations sufficiently high.

Data sets

For the river detection, we tested our method using three DTMs. For some experiments, we only use a smaller section of theses scenes as input data. The sizes and grid size of the data sets are given in Table 1. The synthetic data simulates an ideal tidal channel network and does not include any noise. Moreover, the heights in the channels continuously decrease from the source to the mouth.

The two real data sets are characterized by totally different terrain types. On the one hand, we extracted a tidal channel network from the Wadden Sea. Here, the terrain is nearly horizontal with height differences of 10 m and 1.7 m in the whole scene and the smaller section. In comparison to the synthetic data, the borders of the channels are less clear and the transition to land is smoother, especially for the small channels. On the other hand, we tested our approach on a river network in Vorarlberg, Austria. Here, the terrain is mountainous with height differences of nearl 700 m and 250 m for the whole scene or the smaller test site, respectively. For the synthetic data and the Wadden Sea data, we manually labeled the borders of the channels in the DTM and generated a binary mask of the network. Then, we used a contour-pruned skeleton approach [START_REF] Strzodka | Generalized distance transforms and skeletons in graphics hardware[END_REF] as implemented in Matlab [START_REF] Howe | MATLAB implementation of contour-pruned skeletonization[END_REF] resulting in one pixel wide lines and defined the skeleton as centerlines of the channels. At the head of the streams where standard skeleton algorithms fail we manually digitized the centerlines. For the third data set orthoimages are available which helps to analyze the results. Moreover, reference data of the rivers centerlines are given for this test site. They include some piped river segments which we exclude from the quantitative evaluation.

For the detection of different types of networks in images, we tested our algorithm from the datasets of [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF] that include images of river networks, leaf and retina vessels as well as regular textures as roof tiles (Figures 17 and18).

Apart from the roof tiles, all of these networks have a tree structure. The sizes of the data sets are given in Table 1. Quantitative comparisons with further approaches in the literature are available for these images.

Parameter setting

In our experiments, we set the parameters to values that were determined empirically based on our previous work [START_REF] Schmidt | Network detection in raster data using marked point processes[END_REF]. Unless noted otherwise, they were kept fixed for all tests. In all experiments, we set the initial temperature in equation ( 4) to T t=0 = 10. The number of iterations was set to t = 1 • 10 6 and t = 10 • 10 6 for the synthetic and the real data, respectively. For the whole Austrian test site it was set to t = 40 • 10 6 . The decreasing factor in the cooling scheme was set to d f = 0.99999998. The probabilities for choosing one of the changes and the probabilites of the kernels were set to be equally probable for 455 all kernels (thus, the corresponding ratios in equations 7, 9, 10, 11 and 12 cancel to 1). The further parameters are listed in Table 2. In general, differences in the parameter settings can be explained by the varying appearance of the channels and rivers in the data sets and the size of the test sites.

Evaluation strategy 460

For the evaluation of the results we choose two different methods: a line evaluation for the river networks in the DTM and a pixel-wise comparison for the results of the networks in the images.

The line evaluation is performed based on the method of [START_REF] Wiedemann | Empirical evaluation of automatically extracted road axes[END_REF] by evaluating the edges (neglecting their width) in our graph with reference data of the channel centerline. This is done, first, by defining a buffer surrounding the reference centerline and the extracted edges in the graph, respectively, having a width of b ref and b ext . Then, we tap points on the edges at regular intervals and searched

Correctness CR = result points within b ref sum of all points in the result ,

Quality Q = CR • CP CR + CP -CR • CP . (21)
In order to analyze the geometrical accuracy we calculate the Root Mean Square

Error RMS of the distances d between all result points within b ref and the corresponding reference points by

RM S d = 1 n p np i=1 d 2 i , (23)
where n p is the number of points in the result within b ref .

For the images, ground truth data are available in the form of binary images labeling all pixels within the network. For a pixel-wise comparison, we evaluated all pixels within the segments of the edges (considering their width)

and calculate the completeness and correctness rate as well as the quality of the results. The quality parameters of the result in Figure 8 are given in Table 4. We also compared the logarithmic and the geometrical cooling in the process for this data set. For the logarithmic cooling, the completeness and correctness rates slightly differ, CP = 98.8 %, CR = 94.0 % in comparison to the geometrical cooling, CP = 99.6 %, CR = 88.6 %. The differences can be mainly explained by the difference of the edges at the streams' heads which are slightly beyond the channels borders in the second case. The RMS of the results is 0.66 m and 0.69 m at a grid size of 1 m.

Results

Quality

We also analyzed the importance of each term in our model for the detection of rivers and tidal channels. This is done by successively excluding one of the terms. If we exclude the data term of the high gradient magnitudes at the channel borders (p g = 0), we only detect some network parts (Fig. 9a). In general, the borders of the edges do not correspond to the channel banks. If we do not take into account the homogeneity of the heights perpendicular to the channels (p h = 0), the quality decreases by 5.2 %. By excluding the penalization of overlapping areas (p o = 0), nodes and edges accumulate near the medial axes of the channels where segments with similar widths are sampled at nearly the same positions (Fig. 9b). Thus, the graph does not represent the correct topology of the network. If we do not take into account the connectivity of the graph (p c = 0), the sampling results in several trees which are not connected (Fig. 9c). Due to the data term, the segments still coincide with the borders of the channels rather well. However, we do not end up in a fully connected graph.

If we exclude the term verifying the consistency of the flow directions (p f = 0), the quality decreases by 6.2 %. In general, our results indicate that the chosen model seems to be suitable to describe the extraction problem of river networks.

Reliability

In order to exclude the influence of the generation of random numbers on the result, we repeated the geometrical cooling 100 times. Figure 10 shows the positions of the nodes in all 100 experiments. It can be seen that most of the nodes lie on the centerline of the reference data. In general, the scattering of the nodes positions is low. It is only high at the stream heads, which can be explained by the data term that is near U d = 0 on the nearly horizontal terrain and, thus, does not increase the energy function. Moreover, we observe that the smallest channel is less often detected than the bigger ones. Gradient magnitudes at the borders are much lower than for the bigger channels and are therefore less often accepted in the sampling process. In order to calculate a precision of our result, we determine the standard deviation of all edge positions perpendicular to the channels at regular distance of 1 m for all edges in the 100 results. We only exclude the edge positions at the stream heads. In this way, we obtain a precision of 0.38 m on average that can be attribute to the variability of the random number generation.

We also analyze the energy curves in the experiments. Figure 11 shows the mean energy of all 100 experiments. In general, the performance can be considered as stable. The standard deviation of the energy decreases with increasing iteration. In the last iteration it is 2.6 % of the total energy which can be explained by the smaller channels which are not detected in all results.

The number of nodes and edges varied by about 7 % in the final results of all experiments (see Table 3).

Applicability to different terrain types

We also evaluated our approach on real data: (1) a test site from the Wadden Sea with nearly horizontal terrain and (2) a mountainous test site in Austria. For both data sets we used the hierarchical approach in different resolution steps and, first, detected the bigger streams in the data of lowest resolution and, then, searched for smaller streams in the higher resolution data. Note that we kept most parameters constant in these steps. The evaluation results

for both scenes and two sections of them are given in Table 4. In general, the quality criteria are high. In the Wadden Sea (section) data set the completeness (CP), correctness (CR) and quality (Q) are CP = 72.7 %, CR = 94.4 % and Q = 69.7 %. All bigger channels are completely covered by line segments from the graph (Fig. 12). However, the lower order channels are not completely extracted. Here, the gradient magnitudes at the channels borders are much smaller than for the bigger channels and often show a smooth transition to land, so that the data term does not minimize the energy function in the same way. For most channels, the contours of the channels correspond well to the borders of the graph segments. The RM S of the position is 2.67 px at 1 m pixel size. For the whole Wadden Sea scene, several tidal channels are found in the different resolutions (Figure 13). However, some of the smaller channels are not detected which leads to a relatively low completeness rate of CP = 50.6 %.

The correctness rate is higher, CR = 78.2 %. By using a pixelwise evaluation instead of the line evaluation described in Section 4.4 the completeness and correctness rates increase to CP = 79.2 % and CR = 80.3 %.

For the data set of Austria the first order stream is completely detected in the first hierarchical step as well. It is completely connected apart from one location where the riverbed is very wide and strongly curved (Fig. 14). We also extract many rivers of higher order which are difficult to see in the orthoimages.

The accuracy results for this data set are CP = 76.1 %, CR = 76.1 % and Q = 61.3 %. The RMS of the position is 0.95 px with 1 m pixel size in the DTM. From the river network with a length of 247.3 km in the whole Austrian test site 156.1 km (CP = 64.5 %) are detected (Figure 15). The results show some misclassifications with terrain structures similar to river networks, e.g.

valleys. Moreover, in the urban area on the left site of the scene some of the rivers are not detected.

For all scenes, we observe that the number of extracted rivers increases by increasing the weight β for the data term. Simultaneously, the correctness rate is reduced. These results confirm our model: By integrating prior knowledge about the network with a higher weight (lower value for β) the correctness rate increases. In contrast, the completeness rate increases by increasing the weight β and, thus, giving more weight to the data term.

We also analyze the influence of the prior term of the flow direction U f on the results. By excluding this term, the quality of the results decreases by 0.5 % for the Austrian test site and by 1.1 % for the Wadden Sea scene (see Table 5). The correctness increases by 3.0 % and 1.7 %. The results show that the integration of 3D knowledge -which differentiates our approach from other marked point process such as [START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF] -helps to improve the correctness of the detected network. However, in regard to the correctness and completeness rates, our approach outperforms the flow routing algorithm for both test sites. Especially for the Wadden Sea test site, our correctness rate is much higher (CR f pp = 94.4 %) in comparison to the flow routing algorithm (CR f r = 75.9 %). For the Austrian test site the improvement in the correctness rate is lower, CR f pp = 76.1 % in comparison to CR f r = 75.1 %. However, we achieve a significant higher completeness rate, CP f pp = 76.1 % in comparison to CP f r = 58.9 %. Especially at the stream heads and near to the image borders, our approach achieves better results than the flow routing algorithm, which can be explained by the fact that the catchment area for these pixels is outside the DTM borders.

Transferability to networks in images

We also applied our method to image data showing different types of linenetwork in images. In comparison to the river network detection, the gradients and the homogeneity of gray values (instead of heights) are calculated in the data term. As we did not use a DTM as input data, the term concerning the physical consistency of the network was neglected or set to a low weight for these data sets. Figure 17 shows the results of our method from three different types of images. In general, major parts of the networks are detected by our approach.

620

For the retina network, we end up in one connected tree. In the results of the leaf and retina image, some of the smaller network parts are missing which can be explained by low gradients at their borders. In the tiles image, the network is completely covered by edges of the forest. Contrary to the junction point process proposed by [START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF], our algorithm -constructed for trees -fails to recover 625 correctly the junctions from the tiles image. Here, small gaps occur and the forest is composed of several trees (see Tables 6 and7).

GT [START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF] ours Leaf Retina Tiles In comparison to other methods based on marked point processes or active contours, we achieve comparable results in terms or completeness, correctness and quality (see Table 8). Nevertheless, these methods are composed of many disconnected components as they do not use a graph structure to describe the network (see Figure 18). Here, our approach has advantages as it delivers the topology of the network by using a graph structure. ) with active contour [START_REF] Rochery | Higher order active contours[END_REF] and marked point process [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF][START_REF] Lafarge | Geometric feature extraction by a multimarked point process[END_REF][START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF] approaches from the river image of Figure 18.

Lacoste [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF] Rochery [START_REF] Rochery | Higher order active contours[END_REF] Lafarge [START_REF] Lafarge | Geometric feature extraction by a multimarked point process[END_REF]

Conclusions

In this paper, we presented a stochastic approach based on marked point 635 processes for the automatic extraction of networks in raster data. In our approach, the network is modeled as an undirected, acyclic graph which is iteratively constructed during the optimization process. The approach is evaluated on synthetic data and on two DTM derived from airborne lidar data. Our experiments show that for all data sets, the most relevant tidal channels and rivers are detected, apart from some smaller inlets. In most cases, we end up in one single tree or only a small number of trees. The quality is > 60% in all scenes.

The geometric accuracy is < 1 pixel for the synthetic and the Austrian test site and < 3 pixel for the Wadden Sea data. We also developed a hierarchical detection strategy and searched for channels in different resolutions of the input data. In this way, the parameter space can be significantly reduced. The integration of 3D information in the model improves the correctness rates in all test sites. Our method outperforms a standard flow routing algorithm and achieves significantly higher accuracy results, especially in nearly horizontal terrain. We also proof the transferability of our approach to different types of networks in image data. Here, main parts of the network are detected by our approach. The results are competitive in accuracy in comparison to further approaches based on point processes in the literature.

In the future, we intend to determine parameter by training in order to be less sensitive to the experimental setting. We also plan to investigate on faster optimization schemes by parallelizing the sampling mechanism as proposed by [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF] in case of Markovian energies, or by discretizing the energy formulation into a binary variable optimization problem [START_REF] Pham | Efficient point process inference for large-scale object detection[END_REF]. Such extensions require, however, an entire remodeling of our energy formulation.

Figure 1 :

 1 Figure 1: Forest point processes. (left) A forest point process generates configurations of points where each point l i is associated with the index of its adjacent points (e.g. point index 2 and 5 for l 3 ) and their branch width (e.g. 1 for branch (l 2 , l 3 ), and 2 for (l 4 , l 6 )). (middle) For any configuration, this formalism allows us to connect the points with a unique adjacency graph which is planar and acyclic. (right) The output network is a set of trees represented by connected line-segments, each line-segment having a specified width.

Figure 3 :

 3 Figure 3: The initial example forest consists of three edges {e 1 , e 2 , e 3 } and five nodes {v 1 , v 2 , v 3 , v 4 , v 5 }. Top: In the birth event node v 6 and edge e 4 are added to the forest. They are removed from the forest during the death event. Bottom: If neighboring nodes are not found for the new node v 6 , a second node v 7 and the edge e 4 are added to the graph.

305( 1

 1 ) node translation (2) edge width change (3) edge connectivity change

2 Figure 4 :

 24 Figure 4: In the modification event, the parameters of the forest are changed. Modification (1): We change the coordinates of node v 1 of edge e 1 . This is done by randomly generating a transition vector b in a local neighborhood. The inverse modification is the shift of node v 1 in the opposite direction. Modification (2): We randomly change the width of edge e 1 .Modification (3): The connectivity of node v 5 is changed. Here, we add edge e 4 to the configuration and connect v 5 and v 3 . In the inverse modification, the edge is removed from the forest.

2 Figure 5 :

 25 Figure5: In the merging event, a node is randomly chosen, here v 5 . It is merged with a node in a predefined neighborhood r, here v 3 . Node v 5 is removed from the graph. Conversely, in the split event a node is randomly chosen, here v 3 , and a transition vector b is randomly generated. A new node v 5 is added to the configuration.

Figure 6 :

 6 Figure 6: A typical example of profile of a tidal channel in our test site.

Figure 7 :

 7 Figure 7: For the calculation of the data term, we take into account the DTM gradients of m j = 1 and m j = 2 and the homogeneity of the heights of m j = 3 and m j = 4 of each edge.

Figure 8

 8 Figure 8 shows how the object configuration in the synthetic data set evolves in the course of the sampling procedure, as edges are added to or removed from the configuration or their parameters are changed. After about 760 • 10 3 iterations the tidal channel network is completely covered by segments. The contours of the channels correspond well to the borders of the edge segments and the resultant graph is completely connected.

Figure 8 :

 8 Figure 8: Results for the synthetic data set after 1 • 10 2 , 1 • 10 4 and 1 • 10 6 iterations.

Figure 9 :

 9 Figure 9: Results for the synthetic data set by excluding (left) the data term of the gradient magnitudes at the channel borders (pg = 0), (middle) the penalization of the overlapping (po = 0) and (right) the evaluation of the graph connectivity (pc = 0) where the blue circles indicate the parts of the network which are not connected.

Figure 10 :

 10 Figure 10: Reference centerline of the tidal channel in the synthetic data set (red) and the sampled node (green) for all 100 repetitions of the experiment.

Figure 11 :

 11 Figure 11: The energy is minimized during the sampling process and converges to the global optimum. The figure shows the mean energy of 100 experiments whereby the error bars represent the standard deviation in the respective iterations.

Figure 12 :

 12 Figure 12: Sampling result for the Wadden Sea (section) area with the detected edges (red) of the graph its segments of the first hierarchical step (green) and the second hierarchical step (yellow). The blue line represents the reference centerline.

Figure 13 :

 13 Figure 13: Results for the whole Wadden Sea scene. With increasing brightness of the blue segments, the resolution of the input data increases.

Figure 14 :

 14 Figure 14: Segments of our approach of the first (yellow) and second hierarchical step (red) and the reference network (blue) for the data set of Austria with an orthoimage in the background (β = 0.04).

Figure 15 :

 15 Figure 15: Results (top) and reference (bottom) for the whole Austrian test site (β = 0.04).

Figure 16 :

 16 Figure 16: The extracted edges of our approach (red) compared to the reference network (blue) and the flow routing result, Ca = 25000 m 2 , (yellow) with the input DTM.

Figure 17 :

 17 Figure 17: Comparisons with a junction point process [5] on different network images. Note that our tree-based network representation does not allow the junctions of the Tiles image to be recovered correctly.

Figure 18 :

 18 Figure 18: Comparisons with active contour and marked point process approaches.

Table 1 :

 1 Data sets in our evaluation.

	Data set (DTM)	size of the scene [m 2 ] grid size [m] ∆z[m]
	Synthetic data	170 x 170	1	1.7
	Wadden Sea	8000 x 3500	2, 5, 10, 20	10.0
	Wadden Sea (section)	1001 x 1001	1, 2	1.7
	Austria	10000 x 7500	5, 10	700
	Austria (section)	2500 x 2500	5, 10	250
	Data set (image) size of the scene
		River	364 x 320
		Leaf	651 x 364
		Retina	365 x 378
		Tiles	672 x 523

Table 2 :

 2 Parameter setting in the experiments. If the parameters vary with resolution of the input data, several parameters are listed.

	Parameter Synthetic	Wadden	Wadden Sea	Austria	Austria
		data		Sea	(section)		(section)
	λ	50		500	500	4000	500
	β	0.13		0.15	0.35	0.1, 0.04	0.1, 0.04
	r[px]	16		20, 15	15	5, 6	5, 6
	p h	5		3	3	20, 3	35, 3
	p o	300		200	300	500, 400	500
	p c	100		80, 100	100	250, 150	100
	p f	50		50	50	50	50
	c 1	50	150, 50, 20	50	0	0
	c 2	4		4	4	8	8
		Completeness CP =	ref erence points within b ext sum of all points in the ref erence	,

for corresponding points perpendicular to the centerline within the buffer. The quality of the results is calculated by the completeness and correctness rate which analyze if every object in the scene has been extracted and if the objects are correct extracted, respectively. The combination of both rates is given by the quality and can be used for a ranking of different results. The quality criteria are calculated by

Table 3 :

 3 Variability of the results by repeating the experiment 100 times.

	parameter	mean standard deviation
	number of nodes	37.6	2.6 (6.9 %)
	number of edges	36.1	2.6 (7.2 %)
	total energy U	-4720.4	126.3 (2.6 %)

Table 4 :

 4 Accuracy results (completeness CP, correctness CR, quality Q) by using buffer b for our method (fpp) with varying parameter settings for all test sites.

	Data set	Method	DTM	b	CP	CR	Q	RMS
			[m]	[m]	[%]	[%]	[%]	[px]
	synthetic	fpp (log)	1	3	98.8	94.0	92.9	0.66
	data
		fpp (geom)			99.6	89.0	88.6	0.69
	Wadden Sea	fpp, β = 0.3	1, 2	10	70.1	96.1	68.9	2.75
	(section)
		fpp, β = 0.4			72.7	94.4	69.7	2.67
	Wadden Sea fpp, β = 0.15	2, 5,		50.6	78.2	44.3	0.34
			10, 20
	Austria	fpp, β = 0.04	5, 10	15	76.1	76.1	61.3	0.95
	(section)
		fpp, β = 0.05			82.2	64.4	56.5	0.95
	Austria	fpp, β = 0.04	5, 10		61.4	68.4	47.9	1.04
		fpp, β = 0.05			64.5	63.1	46.9	1.05

Table 5 :

 5 Analysis of the input of the prior term U f concerning the flow direction.

	Data set	Method CP [%] CR [%] Q [%] RMS [px]
	Wadden Sea (section)	with U f	72.7	94.4	69.7	2.67
		without U f	70.3	92.7	66.6	2.56
	Austria (section)	with U f	76.1	76.1	61.3	0.95
		without U f	78.4	73.1	60.8	0.97

Table 6 :

 6 Comparison of the results of our approach (fpp) with a flow routing algorithm (fr) with varying parameter setting for the catchment area Ca.

	Data set	Method	CP [%] CR [%] Q [%] RMS [px]
	Wadden Sea (section)	fpp	72.7	94.4	69.7	2.67
		fr, C a = 7500	72.0	61.6	49.7	4.35
		fr, C a = 12500	68.7	71.3	53.8	4.39
		fr, C a = 25000	60.3	75.9	50.6	4.52
	Austria (section)	fpp	76.1	76.1	61.3	0.95
		fr, C a = 7500	78.1	63.5	53.9	1.01
		fr, C a = 12500	68.6	73.1	54.8	1.01
		fr, C a = 25000	58.9	75.1	49.3	1.07

Table 7 :

 7 Quantitative comparison (correctness CR [%], completeness CP [%] and qualityQ [%]) with a junction point process[START_REF] Chai | Recovering line-networks in images by junction-point processes[END_REF] from the different network images of Figure17.

			Leaf			Retina			Tiles
		CR	CP	Q	CR	CP	Q	CR	CP	Q
	[5]	58.8 70.6 47.3 60.4 54.4 40.1	46.2	64.9	36.9
	Ours 69.5 59.7 47.3	57.2	52.2	37.6	65.1 87.6 59.6

Table 8 :

 8

Quantitative comparison (correctness CR [%], completeness CP [%] and quality Q [%]

As the number of parameters does not change in this event, equation 10 is not a Green ratio in the proper sense, but rather the transition probability of the Metropolis-Hasting algorithm ([START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]).

Acknowledgement

We would like to thank the State Department for Waterway, Coastal and Nature Conservation (NLWKN) of Lower Saxony for providing the Wadden Sea data and the State Office of Surveying and Geoinformation of Vorarlberg, Austria (http://www.vorarlberg.at/lvg), for providing the lidar and the reference data of Vorarlberg.

Comparison with standard flow routing algorithms

We also compared our method with a standard flow routing algorithm and applied the Multiple Flow Direction approach as implemented in GRASS-GIS (r.stream.extract, [START_REF] Jasiewicz | A new GRASS GIS toolkit for Hortonian analysis of drainage networks[END_REF]). This method was also analyzed in detail for a region of the same Austrian data set by [START_REF] Vetter | Modification of high resolution airborne laser scanning DTMs for drainage network delineation, Photogrammetrie, Fernerkundung[END_REF]. In general, flow routing algorithms are computationally efficient; they calculate the flow of the water by exclusively analyzing height differences of neighboring cells in the DTM. That is why the flow routing algorithm outperforms the forest point process regarding the computing time, with computing times of a few seconds in comparison to a few minutes of our approach for the test sites shown in Figure 12 and 14. The quantitative and visual comparison for both test sites are given in Table 4 and Figure 16. We selected the following parameter setting for the flow routing algorithm: we set the minimum segment length to 100 m (similar to [START_REF] Vetter | Modification of high resolution airborne laser scanning DTMs for drainage network delineation, Photogrammetrie, Fernerkundung[END_REF]) and evaluated different settings for the area of the catchment area C a : C a = 7500 m 2 , C a = 12500 m 2 , C a = 25000 m 2 (the latter is similar to [START_REF] Vetter | Modification of high resolution airborne laser scanning DTMs for drainage network delineation, Photogrammetrie, Fernerkundung[END_REF], the other values are chosen to be smaller as our test site is smaller). In general, the completeness rate decreases and the correctness rate increases by increasing the catchment area.

The visual comparison of the results in Figure 16 shows that the flow routing algorithm gives a description of the river centerlines in the form of zigzag lines.

In contrast, the approximation by line segments in the forest point process results in straight lines which represent the river centerlines in a better way and, thus, give a higher geometric accuracy of the results. For the Wadden Sea test site, the geometric accuracy in our results is higher by a factor 1.6 (RM S f pp = 2.67 px, RM S f r = 4.35 px). For the Austrian test site the difference is lower, RM S f pp = 0.95 px vs. RM S f r = 1.01 px. In contrast, the flow routing algorithm delivers better results than the forest point process concerning the connectivity of the network. Whereas all rivers are connected in the results of the flow routing algorithm, our results show few gaps, especially in those parts of the network where a small river flows into a bigger one. In order to overcome this problem, an extension of our model would be required, e.g. by calculating the prior energy term depending on the width of the line segments.

Appendix

In order to define the Jacobian determinant of the transition function in equation ( 4) from the initial configuration X to the new configuration X , the auxiliary variables h and h are introduced. For instance, for the birth event in Figure 3 (top) where a new node (x, y) and the width w for the new edge are added to the configuration, three auxiliary variables are introduced. The further parameters of the initial configuration X remain. Subsequently, the transition can be written by

Then, the Jacobian determinant for the birth event is given by det

For the death event in Figure 3, the transition can be described by

The Jacobian determinant is similar to equation 25.

Analogously, the event of modification and split-and-merge can be described.

In both cases the absolute value of the Jacobian determinant is 1.