217 
the inequality A > |a| must be satisfied. We choose a > 0 so 
that the energy at -1 is less than that at 1. The derivative term en 
sures the smoothness of 0k, producing a narrow interface around 
the boundary dR interpolating between —1 and +1. 
To introduce prior shape information, a nonlocal term is then 
added to give an energy Ep = Eq + ¿’nl (Rochery et al., 2005): 
£nl(0) = 
- 111^ d 2 x d 2 x' d<f>(x) ■ d<f>(x') * , (2) 
where d is the interaction range. This term creates long-range in 
teractions between points of dR (because d(f>R is zero elsewhere) 
using an interaction function, \I>, which decreases as a function of 
the distance between the points. 
In this paper, the interaction function 'I' is taken to be the mod 
ified Bessel function of the second kind of order 0, Kq. This 
choice, as opposed to that used in (Rochery et al., 2005), al 
lows a wide range of stable branch widths. However, it also al 
lows the width of individual branches to fluctuate rapidly, because 
it reduces branch ‘rigidity’. To allow a wide range of branch 
widths while constraining the rate of change of the widths of 
individual branches (without imposing straightness of branches 
as in (Rochery et al., 2005)), and also to favour other important 
properties of directed networks, it is necessary to augment this 
model with an extra field representing a ‘flow’ in the network. 
This leads to the phase field HOAC model for directed networks, 
to be described next. 
2,2 Directed network model 
The phase field HOAC model for directed networks was intro 
duced by (El Ghoul et al., 2009b). To model directed networks, 
the phase field 0 is augmented by a new tangent vector field v 
that loosely speaking represents the ‘flow’ in the network. 1 
The total prior energy Ep(4>,v) is then the sum of a local term 
Eo(<fi,v) and the nonlocal term £nl(0) given by equation (2). 
Eq is given by 
£o(0,v) = J ^ 30 ■ 30 + ^ (d ■ v 
L v 
+ dv : dv + W {(f)., 
(3) 
W((f>, v) is a potential with two degenerate sets of minima, at 
(0, |v|) = (—1,0) and (<f), |v|) = (1,1). These minima define 
the stable phases corresponding to R and R respectively, replac 
ing (f) = ±1 in the undirected model. The potential W is a fourth 
order polynomial in <f> and \v\, constrained to be differentiable: 
W(4>,V) = -Lj (A22 + A21 (f> + A2o)^j— 
+ Aq4 + A03-^- + ^02^- + Aoi(f> ■ 
(4) 
1 To avoid misunderstanding, we stress that v is not intended to rep 
resent the physical flow of, e.g. water, in the network, nor is the model 
intended to model the physical behaviour of the flow. Rather, v is an aux 
iliary variable that acts to favour certain geometric properties of the net 
work. (Since it is not coupled to the image, it is, probabilistically speak 
ing, a ‘hidden variable’.) At the same time, it is not a coincidence that v 
shares many of the properties of physical flows, such as smoothness and 
conservation, nor that the resulting stable configurations resemble physi 
cal flows in the network. 
The second term in equation (3) penalizes the divergence of v. 
This represents a soft version of flow conservation, but in prac 
tice the parameter multiplying this term will be large, so that in 
general the divergence will be small. The third term is a small 
overall smoothing term on v (dv : dv = n (d m v n ) 2 , where 
m,n G {1,2} label the two Euclidean coordinates), since con 
straining the divergence is not sufficient to ensure smoothness. 
Because of the transition from |u| = 1 to |u| = 0 across the 
boundary of the region, the divergence term tends to make v 
parallel to the region boundary, since this results in zero diver 
gence. The smoothness and divergence terms then propagate this 
parallelism to the interior of the branch, with the result that the 
flow tends to be along the branch. This fact, when coupled with 
the constraint on |u| inside the channel, means that width vari 
ations are constrained to be slow along a channel, since total 
flow is directly related to branch width. At the same time, the 
use of = Kq, means that different branches may have very 
different widths. At junctions, the conserved flow along each 
branch favours ‘conservation of width’: the (soft) constraint that 
total incoming flow be approximately equal to total outgoing flow 
translates to the sum of the incoming widths being approximately 
equal to the sum of the outgoing widths. Thus the introduction of 
the new field v can favour network regions with geometric prop 
erties characteristic of directed networks. 
2.2.1 Parameter settings Requiring (—1,0) and (1,1) to be 
extrema of the potential w reduces the number of free parame 
ters of W from seven to four, while requiring these points to be 
minima (i.e. the Hessian at these two points should be positive- 
definite) generates further lower and upper bounds on the remain 
ing parameter values. 
We fix further relations between the parameters by requiring that 
the two minima described above be the only local minima; that W 
be bounded below; and that the potential energy of the network 
region R be greater than of the background R, i.e. w( 1,1) > 
w(—1,0). The resulting potential has a saddle point lying be 
tween the two minima at a point (<f> s ,v s ). This point plays an 
important role: the ‘neutral’ initialization of the gradient descent 
algorithm is given by (0, |u|) = (<f> s ,v s ), the direction of v be 
ing random. In addition, we constrain the parameter /3 in ¿?nl so 
that the part of Ep containing derivatives, i.e. everything except 
W, be positive definite (it is a quadratic form). Since constant 
values of 0 and v produce zero in these derivative terms, which is 
the global minimum value of these terms, and since constant val 
ues of 0 and v equal to those at the global minimum of W, which 
is W(—1,0), produce the global minimum of W, the global min 
imum of Ep is at (0, |u|) = (—1,0). 
The energy Ep can favour different stable geometric configura 
tions depending on the values of the parameters remaining after 
the above constraints have been imposed. Since we are interested 
in modelling networks, we need to choose parameter values that 
favour networks as stable structures. Such values can be found 
using a stability analysis of the model. We assume that network 
branches are long enough and straight enough that their stability 
can be analysed by considering the limit of a long, straight bar, 
whose symmetry facilitates the analysis. We do not detail here 
the stability calculations because they are lengthy: they will be 
reported elsewhere. The stability analysis of a network branch 
places constraints on the parameter values of the model. In all 
the experiments, we use these constraints to fix further parame 
ter values, and to replace others with physical parameters such as 
average branch width.