Ting Peng

Ian H Jermyn

Véronique Prinet

Josiane Zerubia

Extended Phase Field Higher-Order Active Contour Models for Networks Its Application to Road Network Extraction from VHR Satellite Images

Keywords: active contour, phase field, shape prior, parameter analysis, remote sensing, road network extraction

This paper addresses the segmentation from an image of entities that have the form of a 'network', i.e. the region in the image corresponding to the entity is composed of branches joining together at junctions, e.g. road or vascular networks. We present new phase field higher-order active contour (HOAC) prior models for network regions, and apply them to the segmentation of road networks from very high resolution satellite images. This is a hard problem for two reasons. First, the images are complex, with much 'noise' in the road region due to cars, road markings, etc., while the background is very varied, containing many features that are locally similar to roads. Second, network regions are complex to model, because they may have arbitrary topology. In particular, we address a limitation of a previous model in which network branch width was constrained to be similar to maximum network branch radius of curvature, thereby providing a poor model of networks with straight narrow branches or highly curved, wide branches. We solve this problem by introducing first an additional nonlinear nonlocal HOAC term, and then an additional linear nonlocal HOAC term to improve the computational speed. Both terms allow separate control of branch width and branch curvature, and furnish better prolongation for the same width, but the linear term has several advantages: it is more efficient, and it is able to model multiple widths simultaneously. To cope with the difficulty of parameter selection for these models, we perform a stability analysis of a long bar with a given width, and hence show how to choose the parameters of the energy functions. After adding a likelihood energy, we use both models to extract the road network quasi-automatically from pieces of a QuickBird image, and compare the results to other models in the literature. The state-of-theart results obtained demonstrate the superiority of our new models, the importance of strong prior knowledge in general, and of the new terms in particular.

Introduction

The need to segment network-like structures from images arises in a variety of domains. Examples include the segmentation of road and river networks in remote sensing imagery, and of vascular networks in medical imagery. Automatically extracting the network region in the image is a difficult task, because images are usually complex, containing much noise and confounding elements having similar local properties to the entity of interest. For this reason, techniques that include no prior knowledge about the region containing the network cannot succeed. In order to solve this problem, such prior knowledge must be injected somehow, either through the intervention of a user, or by incorporating it into a model. Human users possess very specific prior knowledge about the shape of regions corresponding to networks, and in most applications, this level of knowledge is necessary rather than merely sufficient. Unfortunately, current methods, based on manual extraction, are time and labour intensive. On the other hand, generic prior knowledge alone, for example concerning boundary smoothness, is not enough. The need to include more specific prior knowledge of a class of shapes raises a difficult methodological issue, however. The set of network-like regions is complicated to model, because they may have arbitrary topology. More concretely, it consists of a large (in principle infinite) number of connected components, corresponding to the different possible topologies of a network (number of connected components in the network, number of loops in each connected component), or equivalently to the set of planar graphs (for 2D data). To this is added a geometric superstructure corresponding to an embedding of the graph in the plane, and to its 'fattening' into a region. The construction of a model that favours regions lying in this set as opposed to those outside it is a non-trivial problem. This paper proposes two new models to address this problem, and applies them to the extraction of road networks from very high resolution (VHR) satellite imagery.

The incorporation into models of prior knowledge about a region to be segmented from an image has a long history. The earliest and still most widely used models incorporate local knowledge about the boundary, essentially smoothness: active contours [START_REF] Kass | Snakes: Active contour models[END_REF] are one example, the Ising model [START_REF] Ising | Beitrag zur theorie des ferromagnetismus[END_REF][START_REF] Geman | Stochastic relaxation, Gibbs distribution and the Bayesian restoration of images[END_REF] another. This degree of prior knowledge is almost never sufficient to segment an entity of interest automatically, even in relatively simple images. More recent work has focused on the inclusion of more specific prior knowledge in active contour models [START_REF] Chen | Using prior shapes in geometric active contours in a variational framework[END_REF][START_REF] Cremers | Diffusion snakes: Introducing statistical shape knowledge into the Mumford-Shah functional[END_REF][START_REF] Cremers | Kernel density estimation and intrinsic alignment for shape priors in level set segmentation[END_REF][START_REF] Leventon | Statistical shape influence in geodesic active contours[END_REF][START_REF] Riklin-Raviv | Prior-based segmentation and shape registration in the presence of perspective distortion[END_REF][START_REF] Rousson | Prior knowledge, level set representations & visual grouping[END_REF][START_REF] Srivastava | Statistical shape analysis: Clustering, learning, and testing[END_REF]. This work uses shape priors saying that the region sought must be 'close' to an exemplar region or regions. Although having prior information about the expected shape of the object can significantly increase the robustness of the segmentation algorithm in many applications, this type of model is not appropriate when the region sought has arbitrary topology.

To model families of regions such as networks, [START_REF] Rochery | Higher-order active contours[END_REF] introduced 'higher-order active contours' (HOACs). HOACs incorporate not only local, differential knowledge about the boundary, but also nonlocal, long-range interactions between tuples of contour points. Via such interactions, they favour regions with particular geometric characteristics without constraining the topology via use of a reference region.

For example, the model used in [START_REF] Rochery | Higher-order active contours[END_REF], which uses pairwise interactions, favours, for certain ranges of parameter values, network-like regions composed of branches with roughly parallel borders and a constant width that meet at junctions.

The HOAC energy developed in [START_REF] Rochery | Higher-order active contours[END_REF] suffers from a limitation, however. This is that the interactions between points on the same side of a network branch have the same range and strength as the interactions between points on opposite sides. The effect is that typical maximum curvature κ and branch width W will be related approximately by κ ∼ 1/W . This is particularly serious for certain types of networks, e.g. road networks in cities, for which κ 1/W , and prevents successful extraction.

In this paper, we first construct a new nonlinear nonlocal HOAC prior energy for modelling networks that overcomes this limitation, allowing separate control of branch straightness and width. To optimize computational efficiency, we then propose another new linear nonlocal HOAC prior energy which achieves a similar effect. Moreover, the latter linear energy permits a broader range of widths to be modelled simultaneously, and can even model two disjoint width ranges. We test both models by applying each of them to the problem of road network extraction from VHR images of Beijing. This is an extremely challenging problem due to the amount of 'noise' in the road regions (cars, road markings, shadows, . . . ) and the degree of variation and detail in the nonroad regions. Nevertheless, our new energies permit a quasi-automatic extraction of the road network.

To avoid the complications of expressing regions with arbitrary topology in terms of boundaries and the complexity of the implementation of HOAC terms using standard level-set methods, [START_REF] Rochery | Phase field models and higher-order active contours[END_REF] reformulated HOAC models as equivalent nonlocal phase field models. Phase fields possess many advantages over more traditional methods for region representation and modelling, even in the non-HOAC case, but are particularly advantageous for HOAC energies. It is often convenient to formulate an active contour model in terms of the parametric curve, and then reformulate it as a phase field model for implementation. In this paper, we follow this procedure by considering each new prior energy using first a parametric contour representation, and then a phase field representation.

For the last few decades, road detection from remotely sensed imagery has been extensively studied, due to the variety and importance of the potential applications of an automatic extraction method. A great number of approaches have been proposed [START_REF] Mena | State of the art on automatic road extraction for GIS update: A novel classification[END_REF]. Early work, such as that by [START_REF] Merlet | New prospects in line detection by dynamic programming[END_REF] and [START_REF] Geman | An active testing model for tracking roads from satellite images[END_REF], used path-finding methods suitable for low resolution images, but could only find networks of restricted topology. [START_REF] Péteri | Detection and extraction of road networks from high resolution satellite images[END_REF] took advantage of a topologically correct graph of the network in order to extract roads and junctions using two different types of active contours. [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF] used marked point processes to model road networks, but the models were appropriate for medium resolution images; the method has not been applied to VHR images. [START_REF] Amo | Road extraction from aerial images using a region competition algorithm[END_REF] proposed a region competition based method for providing large-scale GIS information. [START_REF] Hu | Road network extraction and intersection detection from aerial images by tracking road footprints[END_REF] detected roads based on shape classification, and then pruned a road tree using a Bayesian decision process. [START_REF] Bertozzi | Inpainting of binary images using the Cahn-Hilliard equation[END_REF] modified the Cahn-Hilliard equation to achieve fast inpainting of binary imagery. [START_REF] Dobrosotskaya | A wavelet-laplace variational technique for image deconvolution and inpainting[END_REF] developed an extension of this model. Both of them can be used for interpolating simple roads. However, all these methods are restricted to applications on semi-urban areas and using aerial or SAR images. They are not robust enough to be applied to dense urban areas and optical images.

The paper is organized as follows: section 2 recalls HOAC energies and the phase field framework. In section 3, we introduce two new nonlocal HOAC energies. In section 4, we calculate the conditions for which each new total prior model allows stable bars. In section 5, we define the overall model, including a data term. The application of our models to road extraction from VHR images is illustrated in section 6. We conclude in section 7.

2 Higher-Order Active Contours and Phase Fields [START_REF] Rochery | Higher-order active contours[END_REF] proposed a Euclidean-invariant HOAC prior energy for modelling network regions:

E C (R) = λ C L(∂R)+α C A(R)- β C 2 (∂R) 2 dtdt γ(t)• γ(t )Ψ |∆γ(t, t )| d , (1)
where ∂R is the boundary of region R; γ : S 1 → Ω is a map representing ∂R, parameterized by t; Ω ⊂ R 2 is the image domain; dots represent differentiation with respect to t; L is boundary length; A is region area; ∆γ(t, t ) = γ(t) -γ(t ); and d is a constant that controls the range of the interaction. The long range interaction between t and t is modulated by Ψ, the interaction function:

Ψ(x) = 1 2 2 -|x| + 1 π sin(π|x|) if |x| < 2 , 0 else . (2)
It is a smoothly decreasing function from 1 at x = 0 to 0 for x ≥ 2. In equation (1), L(∂R) acts as a regularizer and encourages smoothness of the boundary, while A(R) controls the expansion of the region. The quadratic HOAC term has two effects: it controls the curvature of network branches by trying to align tangent vectors, and it controls branch width by creating a repulsive force. For many reasons [START_REF] Rochery | Phase field models and higher-order active contours[END_REF], the phase field framework provides a more convenient framework for region modelling than do contours. A 'phase field' is a function φ : Ω → R, which defines a region R ∈ Ω via a threshold z: R = ζ z (φ) = {x ∈ Ω : φ(x) > z}. The basic phase field energy term E 0 is

E 0 (φ) = Ω dx 1 2 ∇φ(x) • ∇φ(x) + U (φ(x)) .
(3)

The 'potential' U is given by

U (y) = λ( 1 4 y 4 - 1 2 y 2 ) + α(y - 1 3 y 3 ) , (4)
where λ and α are constants. For λ α > 0, U has two minima, at y = -1 and y = 1, and a maximum at y = α/λ. Define φ R = arg min φ: ζz(φ)=R E 0 (φ). If we ignore the gradient term in equation ( 3), and set z = α/λ, we clearly find that φ R (x) = 1 for x ∈ R and φ R (x) = -1 for x ∈ R = Ω \ R. Adding the gradient term results in a smooth transition from 1 to -1 over an interface region R C around the boundary ∂R.

Note that to a very good approximation ∇φ is non-zero only in R C . It has been shown in [START_REF] Rochery | Phase field models and higher-order active contours[END_REF]) that

E 0 (φ R ) λ C L(∂R) + α C A(R), i.
e. E 0 corresponds to the two linear terms in E C . The third, i.e. the HOAC term in E C , can also be reformulated in terms of an equivalent phase field energy [START_REF] Rochery | Phase field models and higher-order active contours[END_REF]. It becomes

E S (φ) = - β 2 Ω 2 dx dx ∇φ(x) • ∇φ(x ) Ψ |x -x | d . (5)
The sum E 0 + E S is then equivalent to E C in equation ( 1). As explained briefly in section 1, E C (or equivalently E 0 + E S ) suffers from a significant limitation when it comes to modelling networks. Apart from a sign change, the interaction between two points on the same side of a network branch ( γ(t) • γ(t ) > 0) has the same strength and range as the interaction between two points on opposite sides of a network branch ( γ(t) • γ(t ) < 0). The effect is that for a stable network branch, a typical curvature of a branch κ is connected to the width of that branch W approximately by κ ∼ 1/W . In other words, the length/range along which the network branch is expected to be straight is the same as the width of the branch itself. The standard HOAC prior energy E C (or its phase field counterpart E 0 + E S ) thereby provides a poor model of networks with straight narrow branches or highly curved, wide branches. In our application of road extraction in cities, road width gives only an (approximate) upper bound on the radius of curvature of the road: most roads are much straighter than they are wide. For narrow roads, this is particularly problematic, since the road region is relatively unconstrained due to the small road width.

Modelling Networks

To solve these problems, we need to be able to model longer-range or stronger interactions along the road, without changing the interactions across the road. This means that we have to separate the two interaction functions, and hence allow separate control of branch straightness and width, thereby allowing more sophisticated prior knowledge to be included. To achieve this goal, we will construct new nonlocal HOAC prior energies that act in a complementary way to the standard HOAC term. The first new energy term is a nonlinear, nonlocal HOAC energy E N L , which increases the magnitude of the interaction along one side of a network branch. The second is a linear, nonlocal HOAC energy E L , which provides a longer-range interaction along one side of a network branch. Through the stability analysis of this model, we further demonstrate that the linear nonlocal term permits the modelling of two widths simultaneously. These new models have been previously introduced in (Peng et al., 2008b) and (Peng et al., 2008a). Here we present and analyse them in more detail.

Nonlinear Nonlocal HOAC Term

In order to separate the two interactions, the interaction function must depend on the tangent/normal vectors at the pairs of points that are its argument. Although the length scale in the interaction function of equation (1), d, could be made to depend on the inner product between the tangent/normal vectors at the two pixels, it would lead to complicated functional derivatives. Alternatively, we prefer to perform a linear interpolation between two interaction functions. In the contour formulation, our new HOAC prior energy E HO takes the form:

E HO (γ) = - S 1 ×S 1 dsds f + ( γ(s)• γ(s ))Ψ + -f -( γ(s)• γ(s ))Ψ -, (6)
where γ : S 1 → Ω, is an arc length parameterization of the region boundary ∂R; γ(s) is the tangent vector to the boundary at s (thus γ(s) • γ(s ) ∈ [-1, 1]); '+ ' denotes parallel vectors and '-' denotes antiparallel vectors. We define

f + (x), f -(x) : [-1, 1] → [0, 1] by: f + (x) = (1 + x)/2 , (7a) f -(x) = (1 -x)/2 . (7b)
Ψ + and Ψ -are interaction functions similar to that in equation ( 2), but have different range or magnitude. They compete with each other: when γ(s) • γ(s ) ∈ [0, 1], i.e. the two interacting tangent vectors are more parallel, Ψ + is dominant; while when γ(s) • γ(s ) ∈ [-1, 0], i.e. the two interacting tangent vectors are more antiparallel, Ψ -is dominant. Here we decide to adjust only the magnitude of the interaction (although this effectively changes its range also). Thus, we assume that the magnitude of the interaction of parallel vectors is stronger than that of antiparallel vectors, i.e. Ψ + = aΨ -, where a > 1, is a constant. Then, equation ( 6) becomes

E HO (γ) = - 1 2 S 1 ×S 1 ds ds (a -1) + (a + 1) γ(s) • γ(s ) Ψ -. (8)
In order to implement E HO (γ) in the phase field framework, it needs to be reformulated as a function of the phase field φ, instead of the arc length parameterization γ used in equation ( 8). Since the constant length of γ(s) corresponds to the fixed change in φ across the interface, we replace tangent vectors by normal vectors, and then normal vectors by ∇φ. Subsequently, the range of interactions is extended from the region boundary ∂R to the whole of the image domain Ω. Due to the fact that ∇φ(x) is approximately equal to zero everywhere outside the narrow interface R C in Ω, the boundary indicator function

S(φ) = ∇φ(x) • ∇φ(x) ∇φ(x ) • ∇φ(x ) 1 ∀x, x ∈ R C , 0 otherwise , (9)
is inserted into the first term of equation ( 8). Thus we have

E HO (φ) = - 1 2 Ω 2 dxdx (a-1)S(φ)+(a+1) ∇φ(x)•∇φ(x ) Ψ |x -x | d . (10)
When a = 1, this reduces to the standard phase field HOAC term E S (up to a factor of β/2). Therefore, we define our new additional energy term E N L by

E N L (φ) = - β 2 4 Ω 2 dxdx ∇φ(x)•∇φ(x) ∇φ(x )•∇φ(x ) Ψ |x -x | d , (11)
where β 2 > 0 is a constant. (Note that d is the same as that in E S .) Thus in this case, we have a term that is quartic in ∇φ.

The functional derivative of

E N L is δE N L (φ) δφ(x) = β 2 Ω dx ∇ 2 φ(x) Ψ |x -x | d ∇φ(x )• ∇φ(x ) + ∇φ(x ) • ∇φ(x ) ∇φ(x) • ∇Ψ |x -x | d . (12)
Since this functional derivative of E N L contains a term nonlinear in ∇φ, due to S(φ) being O(φ 4 ), we refer to it as the nonlinear nonlocal term.

Whether the two tangent/normal vectors at a pair of interacting points are parallel or antiparallel, the effect of E N L is always to encourage two points inside the range of the interaction to attract each other. Thus E N L weakens the repulsive effect of E S between opposite sides, so that E S along a network branch can be strengthened without changing the width. Consequently, the interaction between pairs of points on the same side of a network branch is stronger than that between pairs of points on opposite sides of a network branch.

Linear Nonlocal HOAC Term

In the previous subsection, we proposed a nonlinear nonlocal HOAC prior term E N L to overcome the limitation of the standard HOAC prior term. The term E N L reinforces the interaction along the bar branch and weakens the interaction across the bar branch, but both interactions still depend on the same interaction function, hence the effects in these two directions cannot be tuned in a completely independent way. In this section, as an alternative to E N L , we construct a new, Euclidean invariant linear nonlocal HOAC prior term E L that works in rather a different way.

As we have seen in equation ( 1), one general class of quadratic HOAC terms can be written as

E HO (γ) = - (∂R) 2 dt dt γ(t) • G C (γ(t), γ(t )) • γ(t ) , (13)
where G C is a map from Ω 2 to 2 × 2 matrices. Imposing Euclidean invariance on this term leads to several possibilities. One is G C (γ(t), γ(t )) = Ψ(|∆γ|/d) δ, where δ is the unit matrix, and ∆γ = γ(t)-γ(t

). Another is G C (γ(t), γ(t )) = Ψ(|∆γ|/d)∆γ∆γ T .
The former leads to the standard HOAC term in E C (equation ( 1)). The latter leads to our new linear nonlocal HOAC prior energy, E L :

E L (γ) = - (∂R) 2 dtdt γ(t)•∆γ(t, t ) × γ(t )•∆γ(t, t ) Ψ |∆γ(t, t )| d 2 , ( 14
)
where we use the same Ψ (equation ( 2)) as in E C , but with a different range d 2 . E L compares each tangent vector to the vector ∆γ(t, t ) joining the two interacting points. When two points have tangent vectors that are both nearly aligned or anti-aligned with ∆γ, the product of the dot products is positive. The energy E L can decrease further by further aligning these tangent vectors with ∆γ and hence with each other. This situation corresponds to two points on the same side of a network branch, as shown in Fig. 1(a). The energy thus favours straight lines, within a range controlled by d 2 . On the other hand, when at least one of the two tangent vectors is nearly orthogonal to ∆γ, the product of the dot products is small. In this configuration, changing the distance between the two points in the argument to Ψ does not change the energy much, and thus the force between two such points is small. This situation corresponds to two points on opposite sides of a network branch, as shown in Fig. 1

(b).

As a result, when E L is added to E C , the width of the network branches is controlled largely by the parameter d of E C , while the distance over which the branch will be straight is controlled largely by d 2 , if d 2 > d. For thin, straight bars, we will indeed fix d 2 > d. The exception to this rule is again shown in Fig. 1(b). From the above, γ(t ) exerts no force on γ(t), but for both γ L (t ) and γ R (t ), the product of the dot products is negative. The energy E L can decrease when the value of Ψ becomes less positive, i.e. both γ L (t ) and γ R (t ) repel γ(t), as shown by the force arrows F L and F R in the figure. The tangential parts of F L and F R cancel, and there is an overall normal repulsion F . If the weight of E L in the model is too large, this repulsion may begin to dominate the bar width.

Making the following change of variables in the linear nonlocal HOAC term E L in equation ( 14):

γ(t) = (cos θ, sin θ), γ(t ) = (cos θ , sin θ ), ∆γ(t, t ) = |∆γ(t, t )|(cos η, sin η) ,
E L can be rewritten (up to a multiplicative factor) as

E L (γ) = - (∂R) 2 dtdt |∆γ(t, t )| 2 cos( θ-θ )+cos θ+ θ Ψ |∆γ(t, t )| d 2 , (15)
where θ = θ -η and θ = θ -η are the angles made by γ and γ with ∆γ. Thus E L is a function not only of θθ , like the standard HOAC quadratic term in equation ( 1), but also of θ + θ . We will continue to use the form of E L given in equation ( 14), however, as this is more convenient for analysis and implementation. when at least one of the two tangent vectors is nearly orthogonal to ∆γ, there is only a very small force between the two points, but contributions from many points can add up to a significant repulsion, as denoted by F

We now reformulate E L (γ) in the phase field framework. We rotate tangent vectors to normal vectors, and replace the latter by ∇φ. Since ∇φ is very small outside R C , the domains of integration can be extended from ∂R to Ω without significantly changing the energy, except for a multiplicative factor. By introducing a weight parameter β 3 , we define the linear nonlocal HOAC phase field term E L (φ) as

E L (φ) = - β 3 2 Ω 2 dxdx ∇φ(x)×(x-x ) ∇φ(x )×(x-x ) Ψ |x -x | d 2 , ( 16
)
where × is the 2D vectorial antisymmetric product. The functional derivative of

E L is δE L (φ) δφ(x) = β 3 Ω dx ∇ • (x -x )(x -x ) T T • ∇φ(x ) Ψ |x -x | d 2 , (17)
where rotates the tangent vectors to the inward normal vectors. Since this functional derivative is linear in φ, we call E L the linear nonlocal term.

Stability Analysis

We now add the new nonlinear term E N L or the new linear term to the standard model By considering the total prior energy of a long, straight bar of a given width, we establish constraints on the para-meters that ensure that a long network branch of the desired width is a stable configuration of the energy functional. An important sideeffect is that some of the (rather abstract) model parameters are effectively replaced by 'physical' quantities, such as bar and interface width, which we can reasonably fix from numerical or application considerations. Note that, in addition, to guarantee the Turing stability of the model, [START_REF] Rochery | Phase field models and higher-order active contours[END_REF] introduced another constraint, which we will not detail here.

E 0 + E S .

Definition of a Bar

Since network branches are locally like straight bars, we can, to a good approximation, analyse the stability of a long straight bar, of length L and width W << L → ∞. This allows us to ignore boundary effects. Denote the width of the interface by w. Such a bar is shown in Fig. 2.

Ideally, we should minimize the prior energy under the constraint that ζ z (φ) = R bar (where R bar denotes the bar region), and then expand around that point to test stability, but this is very difficult. Instead, we take a simple ansatz for φ Rbar , and study its stability in a low-dimensional subspace of function space; the results may be justified a posteriori by numerical experiments. In [START_REF] Rochery | Phase field models and higher-order active contours[END_REF], a similar procedure was followed, the results comparing favourably to those obtained by more sophisticated 'matched asymptotics'. The ansatz is defined as follows. The phase field is given by φ

(x) = 1 for x ∈ R \ R C ; φ(x) = -1 for x ∈ R \ R C , while in R C , φ changes linearly from 1 to -1.
Now we evaluate both prior energies on this ansatz, per unit length of bar. The constraint that the network branch is a local energy extremum requires the first derivatives with respect to w and W to be zero, while the constraint that it is a local minimum requires that the second derivatives must be positive semi-definite. These requirements provide the stability conditions of the model. We evaluate the energy E P,N L on the ansatz, per unit length of bar, which is denoted by e P,N L . Up to an additive constant, it is given by

e P,N L (w, Ŵ ) = 4 3 α Ŵ d + 4 15 λw + 4 w - 16β 2 d w 2 + 4d(β - 2β 2 w 2 ) 2 Ŵ dη η 2 -Ŵ 2 1 -cos(πη) , (18)
where the scaled width Ŵ = W/d. Its first derivatives are set to zero:

∂e P,N L ∂ Ŵ = 4 3 αd -4 Ŵ d(β - 2β 2 w 2 ) 2 Ŵ dη 1 η 2 -Ŵ 2 1 -cos(πη) = 0 , (19a) ∂e P,N L ∂w = 4 15 λ - 4 w 2 + 16β 2 d w 3 2 Ŵ dη 1 η 2 -Ŵ 2 1 -cos(πη) = 0 . ( 19b
)
This is a complicated system of two-variable equations. To simplify the problem, we fix arbitrarily the value of w beforehand (w = 2 ∼ 4), and then solve equation (19a) to get a sub-optimal solution of Ŵ . In this way, there is no need to calculate the second derivative of e P,N L with respect to w, but we still need to ensure that the second derivative of e P,N L with respect to Ŵ is non-negative. It is given by

∂ 2 e P,N L ∂ Ŵ 2 =4d(β - 2β 2 w 2 ) 2 Ŵ dη 1 + π 2 (η 2 -Ŵ 2 ) η 2 -Ŵ 2 cos(πη) -4d(β - 2β 2 w 2 ) ln 2 + 4 -Ŵ 2 Ŵ ≥ 0 . ( 20
)
The above analysis shows that stability of the nonlinear nonlocal HOAC total prior model is related to the scaled control parameters β = β/α and β2 = β 2 /α, and to the scaled width Ŵ = W/d. Therefore, we obtain the parameter constraints:

β - 2 β2 w 2 = 1 3 Ŵ I 1 ( Ŵ ) , (21a)
I 2 ( Ŵ ) ≥ ln 2 + 4 -Ŵ 2 Ŵ , (21b)
where Now we can draw the diagram of the relationship between β, β2 , and Ŵ , which takes the form of a 3D surface. When w is set to 2, an example of such a surface is illustrated in Fig. 3. The left-hand part of the surface represents the extremum as a local maximum, i.e. the second derivative is negative (this is an unstable configuration); the right-hand part of the surface represents the extremum as a local minimum, i.e. the second derivative is positive; and the borderline between these two surfaces, where Ŵ = 0.8798, indicates the situation when two extrema of e P,N L become an inflection point (i.e. ∂e P,N L /∂ Ŵ = 0 and ∂ 2 e P,N L /∂ Ŵ 2 = 0). Clearly, once Ŵ is given, β and β2 have to be chosen in the right-hand part of this surface. For each given width Ŵ , there is a curve that indicates all the possible pairs of scaled parameters β and β2 . Note that, due to the scaling of β and β2 , the dependence on α has been eliminated.

I 1 ( Ŵ ) = 2 Ŵ dη 1 η 2 -Ŵ 2 1 -cos(πη) , (22a)
I 2 ( Ŵ ) = 2 Ŵ dη 1 + π 2 (η 2 -Ŵ 2 ) η 2 -Ŵ 2 cos(πη) . ( 22b

Linear Nonlocal HOAC Term

Similarly, the linear nonlocal HOAC total prior energy per unit length of bar, e P,L , is

e P,L (w, Ŵ ) = 4 3 α Ŵ d + 4 15 λw + 4 w + 4βd 2 Ŵ dη η 2 -Ŵ 2 1 -cos(πη) + 4β 3 d 3 2 d2 Ŵ dη η η 2 -Ŵ 2 2 - η d2 + 1 π sin πη d2 . ( 23
)
As for e P,N L , we find the minimum of e P,L by setting its first derivatives, with respect to w and Ŵ , to zero, while ensuring that its second derivative, with respect to Ŵ , is non-negative. For the derivatives, we find

∂e P,L ∂ Ŵ = 4 3 αd -4β Ŵ d 2 Ŵ dη 1 η 2 -Ŵ 2 1 -cos(πη) - 4β 3 Ŵ d 3 d2 2 d2 Ŵ dη η 2 -Ŵ 2 1 -cos πη d2 , ( 24a
) ∂e P,L ∂w = 4 15 λ - 4 w 2 , ( 24b
)
∂ 2 e P,L ∂ Ŵ 2 = -4βd ln 2 + 4 -Ŵ 2 Ŵ + 4βd 2 Ŵ dη 1 + π 2 (η 2 -Ŵ 2 ) η 2 -Ŵ 2 cos(πη) - 4β 3 d 3 d2 2 d2 Ŵ dη η 2 -2 Ŵ 2 η 2 -Ŵ 2 1 -cos πη d2 . ( 24c
)
It follows that stability is related to the three scaled control parameters β = β/α, β3 = β 3 d 2 /α and d2 = d 2 /d, and also to the scaled width Ŵ = W/d. The parameter constraints are then as follows:

1 -3 β Ŵ I 1 ( Ŵ ) - 3 β3 Ŵ d2 I 3 ( Ŵ ) = 0 , (25a)
βI 2 ( Ŵ ) - β3 d2 I 4 ( Ŵ ) ≥ β ln 2 + 4 -Ŵ 2 Ŵ , (25b)
λ = 15 w 2 , ( 25c
)
where I 1 ( Ŵ ) and I 2 ( Ŵ ) have been defined in equation ( 22); and I 3 ( Ŵ ) and I 4 ( Ŵ ) take the form:

I 3 ( Ŵ ) = 2 d2 Ŵ dη η 2 -Ŵ 2 1 -cos πη d2 , (26a)
I 4 ( Ŵ ) = 2 d2 Ŵ dη η 2 -2 Ŵ 2 η 2 -Ŵ 2 1 -cos πη d2 . ( 26b
)
For w, the constraint is trivial, leading to λ = 15/w 2 . The stable width Ŵ , i.e. the value of Ŵ where a local minimum of e P,L is found, depends on the parameters β, β3 and d2 . The term Ŵ I 1 ( Ŵ ) in equation ( 25a) is a simple curve with one maximum inside the interval [0, 2), so when β3 = 0, for a given β, the number of solutions of this equation is at most two: one corresponds to the local minimum of the energy, the other the local maximum. The additional term Ŵ I 3 ( Ŵ )/ d2 is also a simple curve with one maximum but in the interval [0, 2 d2 ). As a result, when d2 is small, the two curves are mostly overlapped, and equation (25a) possesses the same properties as when β3 = 0: the number of possible solutions is zero or two; while when d2 is large enough so that one curve has non-zero values inside a longer interval than the other, the first constraint changes significantly: the number of possible solutions satisfying ∂E P,L /∂ Ŵ = 0, is zero, two, or four. Consequently, there is a singular point of d2 , where the maximum number of solutions jumps from two to four, and accordingly the stable width(s) that the prior energy can model jump(s) from one to two. However, since equation (25a) involves complicated integrals and the piecewise interaction function, it is non-trivial to obtain analytically the number of its solutions. Therefore, we carry out numerical experiments to find an approximate value for this singular point. We find that such a singular point of d2 (named D2 ) indeed exists, with D2 2.7. If d2 is less than D2 , at most one local minimum can be found. If d2 > D2 , there are three cases, depending on the values of β, β3 and d2 : e P,L has no local minimum; e P,L has one local minimum, with either Ŵ 1 (i.e. W d) or Ŵ d2 (i.e. W d 2 ); or e P,L has two local minima, at Ŵ 1 and Ŵ d2 . The two regimes are illustrated in Figs. 4 and5. For both regimes, the associated stable bar width, as well as the possible behaviours of the energy, are shown in the same figure. Fig. 4 gives an example for d2 = 2 < D2 . In the β3 -β plane, the red region and the green region correspond to the situations when e P,L has no local minimum, and one local minimum respectively. The sole separation curve between the two regions corresponds to the situation when the two extrema of e P,L become an inflection point. This curve demonstrates a qualitative change in the behaviour of e P,L .

The case for d2 > D2 (here d2 = 5.5) is shown in Fig. 5. To make things clearer, we label the three curves by numbers; and we name the widths where the two possible local minima are found ŴMIN1 and ŴMIN2 , and the widths where the two possible local maxima are found ŴMAX1 and ŴMAX2 , by the ascending order of their values, so ŴMAX1 ≤ ŴMIN1 ≤ ŴMAX2 ≤ ŴMIN2 . (In the case of equality, the extrema of the energy merge as one or more inflection points.) We examine the situation at each curve and at each intersection point between two curves:

• Curve 1: ŴMAX1 and ŴMIN1 merge together and become an inflection point. See Fig. 6(c).

• Curve 2: ŴMAX2 and ŴMIN2 merge together and become an inflection point. See Fig. 6(e).

• Curve 3: ŴMIN1 and ŴMAX2 merge together and become an inflection point. See Fig. 6(g).

• Intersection point between curve 1 and curve 2: ŴMAX1 merges with ŴMIN1 and also ŴMAX2 merges with ŴMIN2 . Four extrema of the energy merge as two inflection points. See Fig. 7(a).

• Corner point between curve 1 and curve 3: ŴMAX1 , ŴMIN1 and ŴMAX2 merge together. The energy has a very wide maximum (the inflection point of the three extrema) and a local minimum. See Fig. 7(b). The parameters governing energy stability constitute an M -dimensional parameter space C; the width of the bar is a 1-dimensional state space X; and the energy e P is a smooth function on X parameterized by C. In the linear nonlocal HOAC prior model, the dimension M of the parameter space is 3 (there are three parameters β, β3 and d2 in space C), and so the possible stable states of the 1-dimensional state space should form a swallowtail catastrophe [START_REF] Thom | Structural stability and morphogenesis[END_REF]. The numerical study illustrated by the diagrams in Figs. 4 and 5 confirms this statement.

The variety of possible behaviours is important for applications. As well as being able to model networks with branches of more or less fixed width, but with greater 'stiffness' than provided by the model in [START_REF] Rochery | Higher-order active contours[END_REF], the new linear energy can model two widths at the same time. At certain 'critical points' in parameter space, essentially where pairs of minima merge, it can also potentially model a large range of widths.

Overall Model for Linear Network Extraction

In addition to the prior energy E P (E P,N L or E P,L ), we also need a likelihood energy term linking the region R (which in our case corresponds to the road network) to the data, in our case a VHR optical satellite image. In this section, we will introduce this term, E D , and specify some implementation details.

Overall Energy

The overall energy is the sum of the prior energy E P and the likelihood energy E D :

E(φ; I) = E D (I, φ) + θE P (φ) , (27)
where I : Ω → R is the image, and θ ∈ R + is a constant that balances the contributions of the two terms. E D is given by

E D (I, φ) = - Ω dx φ + (x) ln P + (I(x)) + φ -(x) ln P -(I(x)) . (28)
P ± (I) are models of the histograms of the image intensity, inside (+) and outside (-) the road region. They are both mixtures of Gaussians whose parameters are learned a priori, in a supervised way. The quantities φ ± = (1 ± φ)/2 are, by construction, approximately equal to the characteristic functions of R and R. The likelihood energy is quite weak, in the sense that maximum likelihood classification produces very poor results (see Figs. 9(b) and 11(d)), mainly due to the 'noise' in the road region and the great variations in the background. No image model that only takes into account the radiometry of independent pixels can do much better than this, which is why a powerful geometric prior model is needed.

Optimization and Parameter Settings

To minimize E, we perform gradient descent with the neutral initialization: the initial value of φ is set equal to the threshold z = α/λ everywhere in Ω [START_REF] Rochery | Phase field models and higher-order active contours[END_REF]. The algorithm is thus quasi-automatic. The functional derivatives of the HOAC terms δE S /δφ, δE N L /δφ and δE L /δφ involve convolutions: they are calculated in the Fourier domain, as are all derivatives. The evolution equation for E with

E N L is ∂φ(x) ∂t = 1 2 ln P + P - + θ ∇ 2 φ(x) -λ(φ 3 (x) -φ(x)) -α(1 -φ 2 (x)) + βF -1 k 2 d Ψ(kd) φ(k) -β 2 ∇ 2 φ(x)F -1 d Ψ(kd)F ∇φ(x) • ∇φ(x) -β 2 ∇φ(x) • ∇ F -1 d Ψ(kd)F {∇φ(x) • ∇φ(x)} , (29)
where F and F -1 denote the Fourier and the inverse Fourier transform respectively, and a hat ˆindicates the Fourier transform of a function. The evolution equation for E with E L involves replacing the last two lines of equation ( 29) by

β 3 F -1 k • F Ψ |x| d 2 x x T T • k φ(k) . (30)
The time evolution of φ uses the forward Euler method. The parameters of the two prior energies E P,N L and E P,L are constrained by the stability analysis of section 4.

Experimental Results

As input data I, we use a number of images, with average size 1200 × 1200 pixels, extracted from a QuickBird optical panchromatic image of Beijing. The scenes are characteristic of dense urban regions. Our aim is to extract, completely and accurately, the road network from an image. We also analyse the effect of the different terms in our energies. In order to evaluate the performance of our new nonlinear and linear models, we compare them quantitatively to ground truth and to other methods from the literature. Ground truth was created by slightly correcting by hand an available GIS map from a few years earlier of the road network in the zone shown in the image. Note that this ground truth defines what we mean by the phrase 'road network'.

Nonlinear Nonlocal Overall Model

In this subsection, we demonstrate the behaviour of our new model E = θ(E 0 + E S + E N L ) + E D via experiments at reduced resolutions. Fig. 8(a) shows one of the input images at 1/4 resolution. The parameters (θ, α, λ, β, β 2 , d) are (100, 0.12, 3.8, 0.0375, 0.0338, 4). The results obtained using the energy with the new nonlinear nonlocal term E N L at iterations 1, 1, 500 and 27, 000 are illustrated in . The result obtained using the standard HOAC model without E N L (i.e. β 2 = 0) is shown in Fig. 9(a). We see that adding E N L enables the recovery of the main and secondary road network, whereas the model without E N L misses a secondary road. In order to illustrate the effects of other terms in the model, we compute results using maximum likelihood estimation (MLE, i.e. θ = 0) and a standard, nonhigher-order active contour (i.e. β = β 2 = 0) (see Figs. 9(b) and 9(c)). The ground truth is presented in Fig. 9(d). The MLE result shows that local image information alone is not sufficient to distinguish the roads from the background, while the standard active contour result shows the importance of the geometric knowledge introduced by HOACs. Quantitative evaluations based on standard criteria [START_REF] Heipke | Evaluation of automatic road extraction[END_REF] are Table 1. Quantitative evaluation criteria for the different methods as tested on Fig. 8(a) at 1/4 resolution (T = True, F = False, P = Positive, N = Negative). The completeness is the percentage of ground truth road network that is extracted; the correctness is the percentage of extracted road network that is correct; and the quality is the most important measure of the 'goodness' of the result, because it takes into account the completeness and the correctness [START_REF] Wang | Extraction of urban road network using Quickbird pan-sharpened multispectral and panchromatic imagery by performing edge-aided post-classification[END_REF] and [START_REF] Yu | A novel two-steps strategy for automatic GISimage registration[END_REF] shown in Table 1. Fig. 10 presents more results at reduced resolutions. The first column shows the input image data, which is either at 1/4 or 1/2 resolution. The two columns on the right show the corresponding results obtained with and without the new nonlinear, nonlocal term E N L . The importance of E N L is clear: it facilitates greatly the retrieval of secondary roads. However the nonlinear nature of the energy means that it is prohibitively expensive computationally at full resolution. In the next subsection, we show results at both reduced and full resolution using the linear nonlocal model.

X X X X X X X X X X Method

Linear Nonlocal Overall Model

As in the previous subsection, we study primarily the extraction of a network consisting of roads of roughly the same width, but in section 6.2.2, we briefly consider the extraction of networks containing roads of two very different widths. In the former case, we choose the parameters so that e P,L has one local minimum. The resulting model can extract roads whose widths' are close to the minimizing value. In the latter case, we choose the parameters so that e P,L has two local minima. Again a small range of widths around each minimum is possible.

Extraction of Roads of Similar Width

We apply our new linear model E = θ(E 0 + E S + E L ) + E D to both full-resolution and reduced resolution images. We fix the parameters as described in subsection 4.3. For all experiments, the parameters (θ, α, λ, β, β 3 , d, d 2 ) are (200, 0.15, 4, 0.02, 2 × 10 -4 , 4, 12) and (200, 0.15, 4, 0.02, 1.25 × 10 -5 , 16, 48) at 1/4 and full resolution respectively. Note that apart from the obvious scaling of d and d 2 , and a change in β 3 , the other parameters are the same for the two resolutions.

The results obtained using the new linear model E (with E L ), at 1/4 resolution and at full resolution, are shown in Figs. 11(b) and 11(c). The complete road network is retrieved successfully at both resolutions. Although the segmentation at 1/4 resolution appears geometrically smoother, the extraction result is actually more accurate at full resolution. Accuracy at 1/4 resolution is limited both directly, by the low resolution of the phase field, and indirectly, because each scaling coefficient in the data at level 2 is Table 2. Quantitative criteria for the experiments shown in Fig. 11(a) at full resolution (except first row) (T = True, F = False, P = Positive, N = Negative). See Table 1 for an explanation of completeness, correctness and quality. To evaluate the performance of the linear model, we now compare our result to other methods at full resolution: with MLE (i.e. θ = 0); with a standard, non-higherorder active contour (i.e. β = β 3 = 0); and with the standard model without E L (i.e. β 3 = 0). The results are shown in . Again MLE discriminates poorly between the roads and the background, while the models without E L and/or E S are not able to recover the complete road network (although that with E S does better than the standard active contour, which has only local prior knowledge). In addition, we apply two other methods, as described in [START_REF] Wang | Extraction of urban road network using Quickbird pan-sharpened multispectral and panchromatic imagery by performing edge-aided post-classification[END_REF] and [START_REF] Yu | A novel two-steps strategy for automatic GISimage registration[END_REF], and compare them to ours (see Fig. 12). [START_REF] Wang | Extraction of urban road network using Quickbird pan-sharpened multispectral and panchromatic imagery by performing edge-aided post-classification[END_REF] presented a classification, tracking, and morphology algorithm to extract urban road networks from QuickBird images; [START_REF] Yu | A novel two-steps strategy for automatic GISimage registration[END_REF] proposed a fast but rough segmentation technique based on 'straight line density' to extract urban road networks. However, without much prior geometric knowledge, both of them extract many incorrect areas that happen to have statistical properties similar to roads. Moreover, the accuracy of the delineation of the road boundary is poor. Some quantitative evaluations based on standard criteria [START_REF] Heipke | Evaluation of automatic road extraction[END_REF], are shown in Table 2. The 'quality' is the most important measure because it considers both completeness and correctness. Our complete model outperforms all others.

X

X X X X X X X X

Extraction of Roads of Different Widths

Images containing roads of different widths are processed after choosing parameter values for which e P,L has two local minima. Fig. 15(a) shows an input image containing two roads: their widths are approximately 20 pixels and 80 pixels. The results obtained using the new linear model (with E L ) and the standard model (without E L ), at 1/4 resolution, are illustrated in Figs. 15(b) and 15(c) respectively. The parameter values (θ, α, λ, β, β 3 , d, d 2 ) used in this experiment are (25, 0.15, 5, 0.02, 1.228 × 10 -4 , 4, 22). The estimated stable widths for these parameter values are 5.28 and 20.68, corresponding to the road widths at 1/4 resolution, i.e. 5 pixels and 20 pixels. This comparison shows clearly that adding E L enables the detection of roads with both widths, while the standard model without E L finds only an incomplete network.

Discussion

Computational cost: The computational costs of the two proposed algorithms depend on two factors: the number of iterations and the number of operations at each iteration. The number of iterations depends greatly on the complexity of the image, and on the values of the parameters, even when the final result does not. On the other hand, the number of operations at each iteration is related to the functional derivative of the energy. As already mentioned, we compute the functional derivative in the Fourier domain. At each iterative step, for the nonlinear model, we need to compute two forward Fourier transforms and two inverse Fourier transforms; while for the linear model, only one forward and one inverse transform are required, the rest of the operations being the same. The computation time for the result in Fig. 8, obtained with E N L , is around 80 minutes, while that for the result in Fig. 11(b), obtained with E L , is around 60 minutes. These computation times are long of course, but to our knowledge no other method even approaches the quality of the results obtained.

Parameter choice: We note that in practice, the results obtained are not sensitive to the precise choice of parameter values, provided they lie in the correct subset of the β -β3 -d2 diagram.

Conclusion

Using the phase field HOAC framework, we have presented, in this paper, two novel nonlocal HOAC energies for modelling network-shaped regions. Both of them enable the generation of longer, straighter branches with better prolongation, but by different means. The first, nonlinear energy causes pairs of points inside the range of the interaction to attract each other. In conjunction with the standard HOAC geometric term, it allows the interaction between points on the same side of a network branch to be stronger than the interaction between points on opposite sides of a network branch. The second linear energy, which includes a longer-range interaction, controls the pro-longation of the network branch, but has little effect on the branch width. The standard HOAC geometric term still controls the branch width with a (relatively) short-range interaction. Based on a stability analysis of a bar with a desired width, we established constraints on the parameters of the energy function. We explored the possible behaviours of each of the two resulting prior energies as a function of the parameter settings. We showed that as well as separating the interactions between points on the same and opposite sides of a network branch, the new linear model permits the modelling of two widths simultaneously. The analysis also fixes some of the model parameters in terms of network width(s). Moreover, the linear nonlocal term is more efficient from a computational point of view, and can thus be applied to images at full resolution. As a consequence, narrower network branches can be extracted than is possible at reduced resolution, and in general, extraction accuracy is improved. Experiments also demonstrate the superiority of our new models to others in the literature. Our current work is focused on constructing a prior energy E P that has a very flat local minimum in a wide range, instead of two sharp local minima. This might be a better solution for the extraction of roads with multiple widths. We note finally that also though in this paper we have focused on reducing curvature, one could use the flexibility of the new models to allow more curvature than the model of [START_REF] Rochery | Higher-order active contours[END_REF], while still allowing a wide range of branch widths.

Figure 1 .

 1 Figure 1. The effects of E L in equation (14) (blue arrow: vector; cyan arrow: interaction force; black dot: interacting point). (a): when two tangent vectors are nearly aligned or anti-aligned with ∆γ, the energy E L favours their alignment; (b):when at least one of the two tangent vectors is nearly orthogonal to ∆γ, there is only a very small force between the two points, but contributions from many points can add up to a significant repulsion, as denoted by F

Figure 2 .

 2 Figure 2. A bar of length L → ∞ and width W

)Figure 3 .

 3 Figure 3. Diagram of the scaled control parameters β and β2 , and the scaled width Ŵ , when w = 2. The left part of the surface represents the parameter space for which e P,N L has one local maximum; in the right part of the surface, e P,N L has one local minimum; and the borderline between them, where Ŵ = 0.8798, represents the situation when two extrema become an inflection point

Figure 4 .Figure 5 .

 45 Figure 4. d2 = 2 < D2 . (a): different regions in the β3 -β plane. e P,L has either no local minimum (red) or one local minimum (green). (b): the associated stable bar width Ŵ * . (c): e P,L with no local minimum ( β = 0.05, β3 = 0.04). (d): e P,L with one local minimum ( β = 0.2, β3 = 0.1)

  Figure 6. Evolution of e P,L as the parameter setting is moved from left to right along the dotted path P shown in Fig. 5 ( d2 = 5.5). Red dot: local minimum

Figure 8 .Figure 9 .

 89 Figure 8. Data, and experiment at 1/4 resolution. (a): image data (size = 350 × 350, road width = 3 ∼ 5 pixels). (b)-(d): results obtained using the new nonlinear model (with E N L ) at iterations 1, 1, 500 and 27, 000

Figure 10 .

 10 Figure 10. More experiments at reduced resolutions. First column: input images, first row: 1/2 resolution, size = 400 × 440, road width = 2 ∼ 4 pixels; second row: 1/4 resolution, size = 300×300, road width = 3 ∼ 5 pixels; third row: 1/4 resolution, size = 300 × 400, road width = 3 ∼ 6 pixels; last row: 1/4 resolution, size = 512 × 512, road width = 3 ∼ 15 pixels. Two rightmost columns: corresponding results obtained using the new nonlinear model (with E N L ) and the standard model (without E N L )

Figure 11 .Figure 12 .

 1112 Figure 11. Experiments using the new linear model (with E L ), and analysis of the effects of the different terms in the energy. (a): image data (size: 1400 × 1400); (b): result obtained using the new linear model (with E L ), at 1/4 resolution; (c): result obtained using the new linear model (with E L ), at full resolution; and results obtained at full resolution using (d): MLE; (e): the model with β = β 3 = 0 (equivalent to a standard active contour); (f): the standard model without E L (i.e. β 3 = 0)

Figure 13 .

 13 Figure 13. More results using the new linear model (with E L ) on pieces of a QuickBird image at full resolution. Image size: top left: 800×880; top right: 1200×1200; bottom left: 1200 × 1600; bottom right: 1400 × 1400

Figure 15 .

 15 Figure 15. Extraction of a road network containing two different widths, at 1/4 resolution. (a)-(c): image data; results using: the new linear model E; the standard model (β 3 = 0)

  The sum of the three energies E P,N L = E 0 + E S + E N L constitutes the nonlinear nonlocal HOAC total prior model and E P,L = E 0 + E S + E L constitutes the linear nonlocal HOAC total prior model. For both models, there are a number of parameters to tune, (α, λ, β, β 2 , d) for the former and (α, λ, β, β 3 , d, d 2 ) for the latter. Unfortunately, not all parameter values allow stable network structures. The behaviour of the prior energy depends on the parameter settings, and can vary significantly. If we wish to model networks with this energy, it is therefore very important to analyse the stability of a network structure and to deduce the resulting constraints on the parameters.

Acknowledgements

This work was partially supported by European Union Network of Excellence MUS-CLE (FP6-507752) and by INRIA Associated Team 'SHAPES'. This work was also supported by the 863 program of the Ministry of Science and Technology of China. The work of the first author was supported by an MAE/Thales Alenia Space/LIAMA grant during her PhD.