A higher-order active contour model of a ‘gas of circles ’ and

We present a model of a 'gas of circles': regions in the image domain composed of a unknown number of circles of approximately the same radius. The model has applications to medical, biological, nanotechnological, and remote sensing imaging. The model is constructed using higher-order active contours (HOACs) in order to include non-trivial prior knowledge about region shape without constraining topology. The main theoretical contribution is an analysis of the local minima of the HOAC energy that allows us to guarantee stable circles, fix one of the model parameters, and constrain the rest. We apply the model to tree crown extraction from aerial images of plantations. Numerical experiments both confirm the theoretical analysis and show the empirical importance of the prior shape information.


Introduction
Forestry is a domain in which image processing and computer vision techniques can have a significant impact. Resource management and conservation require information about the current state of a forest or plantation. Much of this information can be summarized in statistics related to the size and placement of individual tree crowns (e.g. mean crown area and diameter, density of the trees). Currently, this information is gathered using expensive field surveys and time-consuming semiautomatic procedures, with the result that partial information from a number of chosen sites frequently has to be extrapolated. An image processing method capable of automatically extracting tree crowns from high resolution aerial or satellite images and computing statistics based on the results would greatly aid this domain.
The tree crown extraction problem can be viewed as a special case of a general image understanding problem: the identification of the region R in the image domain Ω corresponding to some entity or entities in the scene. In order to solve this problem in any particular case, we have to construct, even if only implicitly, a probability distribution on the space of regions P(R|I, K). This distribution depends on the current image data I and on any prior knowledge K we may have about the region or about its relation to the image data, as encoded in the likelihood P(I|R, K) and the prior P(R|K) appearing in the Bayes' decomposition of P(R|I, K) (or equivalently in their energies − ln P(I|R, K) and − ln P(R|K)). This probability distribution can then be used to make estimates of the region we are looking for.
In the automatic solution of realistic problems, the prior knowledge K, and in particular prior knowledge about the 'shape' of the region, as described by P(R|K), is critical. The tree crown extraction problem provides a good example: particularly in plantations, R takes the form of a collection of approximately circular connected components of similar size. There is thus a great deal of prior knowledge about the region sought. The question is then how to incorporate such prior knowledge into a model for R. If the model does not include enough prior knowledge, it will be necessary for the user to provide it.
The simplest prior information concerns the smoothness of the region boundary.
For example, the Ising model and many active contour models [1,2,3] use a combination of region boundary length and region area as their prior energies, but curvature can be used too [1]. Such models are integrals over the region boundary of a function of various derivatives of the boundary. In consequence, they capture local differential geometric information, corresponding to local interactions between boundary points, but can say nothing more global about the shape of the region.
To go further, one must introduce longer range interactions. There are two principal ways to do this: one is to introduce hidden variables, given which the original variables of interest are (more or less) independent. Marginalizing over the hidden variables then introduces interactions between the original variables. Another is to include explicit long-range interactions between the original variables.
The first approach has been much investigated, in the form of template shapes and their deformations. Here a probability distribution or an energy is defined based on a distance measure of some kind between regions. One region, the template, is fixed, while the other is the variable R. Template regions may be learned from examples [4,5,6,7,8,9] or fixed by hand [10]; similarly the distance function maybe based, for example, on the learned covariance of a Gaussian distribution [5,6,7,8,9], or chosen a priori [4,10,11]. The most sophisticated methods use the kernel trick to define the distance as a pullback from a high-dimensional space, thereby allowing more complex behaviours [12]. Multiple templates may also be used, corresponding to a mixture model [12,13].
These methods assign high probability to regions 'close' to certain points in the space of regions. The set of regions with high probability is thus in some sense bounded. As such, it is difficult to construct models of this type that favour regions for which the topology, and in particular the number of connected components, is unknown a priori, because the set of regions in this case is unbounded, and cannot be described as variations around one or more templates. There are many problems, however, for which the topology is unknown a priori, for example, the extraction of networks, or the extraction of an unknown number of objects of a particular type from astronomical, biological, medical, or remote sensing images. For this type of prior knowledge, a different type of model is needed. Higher-order active contours (HOACs) are one such category of models.
HOACs [14] take the second approach mentioned above. They introduce explicit long-range interactions between region boundary points via energies that contain multiple integrals over the boundary, thus avoiding the use of template shapes.
HOAC energies can be made intrinsically Euclidean invariant, and, as required by the above analysis, incorporate sophisticated prior information about region shape without necessarily constraining region topology. As with other methods incorporating significant prior knowledge, it is not necessary to introduce extra knowledge via an initialization close to the target region: a generic initialization suffices, thus rendering the method quasi-automatic. [14] applied the method to road extraction from satellite and aerial images using a prior which favours network-like objects.
In this paper, we describe a HOAC model of a 'gas of circles': the model favours regions composed of an a priori unknown number of circles of a certain radius. For such a model to work, the circles must be stable to small perturbations of their boundaries, i.e. they must be local minima of the HOAC energy, for otherwise a circle would tend to 'decay' into other shapes. The main theoretical contribution of this paper is an analysis of the stability of local minima of the HOAC energy that allows us to ensure that circles of a given radius are stable. In addition, it allows us to fix one of the model parameters in terms of the others, and to constrain the rest. The type of calculation has wide applicability to other active contour models and to other shapes. For example, it shows that no stable circle is possible using a classical active contour model containing only boundary length and interior area terms. The calculation proceeds by performing a functional Taylor expansion of the HOAC energy around a circle (or more generally, any shape), and then demanding that the first order term be zero for all perturbations, and that the second order term be positive semi-definite. Gradient descent experiments using the HOAC energy, with parameters fixed using the stability calculations, produce stable circles of the expected radii, thereby demonstrating empirically the coherence between the stability calculations and the numerical computations used in practice to minimize the energy.
The model has many potential applications, to medical, biological, physical, and remote sensing imagery in which the entities to be identified are circular. We choose to apply it to the problem of extracting tree crowns from aerial imagery, using the 'gas of circles' model as a prior energy, and an appropriate likelihood. We will see that the extra prior knowledge included in the 'gas of circles' model permits the separation of trees that cannot be separated by simpler methods, such as maximum likelihood or classical active contours. We focus on images of plantations and orchards, for which the model is well adapted. The case of general forests is much harder, and will be left for future work.
In the next section, we present a brief introduction to HOACs. In section 3, we describe the 'gas of circles' HOAC model, the stability analysis, and the results of geometric experiments. In section 4, we apply the new model to tree crown extraction. We describe a likelihood energy for trees, and then present experimental results on synthetic data and on aerial images. We conclude in section 5, and discuss some open issues with the model.

Higher order active contours
HOAC models, like all active contour models, represent a region R by its boundary, ∂R, a closed 1-chain γ in the image domain Ω ( [15] is a useful reference for the following discussion). Although region boundaries correspond to a special subset of closed 1-chains known as domains of integration, active contour energies themselves are defined for general 1-chains. It is convenient to use this more general context to distinguish HOAC energies from classical active contours, because it allows for notions of linearity to be used to characterize the complexity of energy functionals.
Using this representation, HOAC energies can be defined as follows [14]. Let γ be a 1-chain in Ω, and dom γ be its domain. Then γ n : (dom γ) n → Ω n is an n-chain in Ω n . We define a class of (n − p)-forms on Ω n that are 1-forms with respect to (n − p) factors and 0-forms with respect to the remaining p factors (by symmetry, it does not matter which p factors). These forms can be pulled back to (dom γ) n by γ n . The Hodge duals of the p 0-form factors with respect to the induced metric on dom γ can then be taken independently on each such factor, thus converting them to 1-forms, and rendering the whole form an n-form on (dom γ) n . This n-form can then be integrated on (dom γ) n .
In the (n, p) = (n, 0) cases, we are simply integrating a general n-form on the image of γ n in Ω n , thus defining a linear functional on the space of n-chains in Ω n , and hence an n th -order monomial on the space of 1-chains in Ω. Taking arbitrary linear combinations of such monomials then gives the space of polynomial functionals on the space of 1-chains. By analogy we refer to the general (n, p) cases as 'generalized n th -order monomials' on the space of 1-chains in Ω, and to arbitrary linear combinations of the latter as 'generalized polynomial functionals' on the space of 1-chains in Ω. HOAC energies are generalized polynomial functionals.
Standard active contour energies are generalized linear functionals on 1-chains in this sense, hence the term 'higher-order'.
The (1, 1) case is simply the boundary length in some metric. The (1, 0) case gives the region area in some metric. An interesting application of the (2, 2) case to topology preservation is described by [16]. We specialize to the (2, 0) case. Let F be an 2-form on Ω n . Using the antisymmetry of F together with the symmetry of γ 2 , we can write the energy functional in this case as and τ =γ is the tangent vector to γ.
By imposing Euclidean invariance on this term, and adding linear terms, [14] defined the following higher-order active contour prior: where L is the boundary length functional, A is the interior area functional and R(t, t ′ ) = |γ(t) − γ(t ′ )| is the Euclidean distance between γ(t) and γ(t ′ ). [14] used the following interaction function Φ: In this paper, we use this same interaction function with d = ǫ, but other monotonically decreasing functions lead to qualitatively similar results.

The 'gas of circles' model
For certain ranges of the parameters involved, the energy (2.2) favours regions in the form of networks, consisting of long narrow arms with approximately parallel sides, joined together at junctions, as described by [14]. It thus provides a good prior for network extraction from images. This behaviour does not persist for all parameter values, however, and we will exploit this parameter dependence to create a model for a 'gas of circles', an energy that favours regions composed of an a priori unknown number of circles of a certain radius.
For this to work, a circle of the given radius must be stable, that is, it must be a local minimum of the energy. In section 3.1, we show that stable circles are indeed possible provided certain constraints are placed on the parameters. More specifi-cally, we expand the energy E g in a functional Taylor series to second order around a circle of radius r 0 . The constraint that the circle be an energy extremum then requires that the first order term be zero, while the constraint that it be a minimum requires that the operator in the second order term be positive semi-definite. These requirements constrain the parameter values. In subsection 3.2, we present numerical experiments using E g that confirm the results of this analysis.

Stability analysis
We denote a member of the equivalence class of maps representing the 1-chain defining the circle by γ 0 , and a small perturbation by δγ. To second order, where ·|· is a metric on the space of 1-chains.
Since γ 0 represents a circle, it is easiest to express it in terms of polar coordinates r, θ on Ω. For a suitable choice of coordinate on S 1 , a circle of radius r 0 centred on the origin is then given by γ 0 (t) = (r 0 (t), θ 0 (t)), where r 0 (t) = r 0 , θ(t) = t, and t ∈ [−π, π). We are interested in the behaviour of small perturbations δγ = (δr, δθ). Because the energy E g is defined on 1-chains, tangential changes in γ do not affect its value. We can therefore set δθ = 0, and concentrate on δr.
On the circle, using the arc length parameterization t, the integrands of the different terms in E g are functions of t − t ′ only; they are invariant to translations around the circle. In consequence, the second derivative δ 2 E g /δγ(t)δγ(t ′ ) is also translation invariant, and this implies that it can be diagonalized in the Fourier basis of the tangent space at γ 0 . It is thus easiest to perform the calculation by expressing δr in terms of this basis: δr(t) = k a k e ir 0 kt , where k ∈ {m/r 0 : m ∈ Z}. Below, we simply state the resulting expansions to second order in the a k for the three terms appearing in equation (2.2). Details can be found in appendix A.
The boundary length and interior area of the region are given to second order by Note that there are no stable solutions using these terms alone. For the circle to be an extremum, we require λ C 2π + α C 2πr 0 = 0, which tells us that α C = −λ C /r 0 .
The criterion for a minimum is, for each k, λ C r 0 k 2 + α C ≥ 0. We must have λ C > 0 for stability at high frequencies. Substituting for α C , the condition becomes The quadratic term can be expressed to second order as π −π dt dt ′ G(t, t ′ ) = 2π π −π dp F 00 (p) + 4πa 0 π −π dp F 10 (p) where Note that there are no off-diagonal terms linking a k and a k ′ for k = k ′ : the Fourier basis diagonalizes the second order term.

Parameter constraints
Note that a circle of any radius is always an extremum for non-zero frequency perturbations (a k for k = 0), as these Fourier coefficients do not appear in the first order term (this is also a consequence of invariance to translations around the circle). The condition that a circle be an extremum for a 0 as well (e 1 = 0) gives rise to a relation between the parameters: where we have introducedr 0 to indicate the radius at which there is an extremum, to distinguish it from r 0 , the radius of the circle about which we are calculating the expansion (3.1). The left hand side of figure 1 shows a typical plot of the energy e 0 of a circle versus its radius r 0 , with the β C parameter fixed using the equation (3.6) with λ C = 1.0, α = 0.8, andr 0 = 1.0. The energy has a minimum at r 0 =r 0 as desired. The relationship betweenr 0 and β C is not quite as straightforward as it might seem though. As can be seen, the energy also has a maximum at some radius.
It is not a priori clear whether it will be the maximum or the minimum that appears Right: the second derivative of E g , e 2 , plotted againstr 0 k for the same parameter values.
The function is non-negative for all frequencies.
atr 0 . If we graph the positions of the extrema of the energy of a circle against β C for fixed α C , we find a curve qualitatively similar to that shown in figure 2 (this is an example of a fold catastrophe). The solid curve represents the minimum, the dashed the maximum. Note that there is indeed a unique β C for a given choice ofr 0 .
Denote the point at the bottom of the curve by (β the extrema merge and for β C < β (0) C , there are no extrema: the energy curve is monotonic because the quadratic term is not strong enough to overcome the shrinking effect of the length and area terms. Note also that the minimum cannot move below r 0 = r (0) 0 . This behaviour is easily understood qualitatively in terms of the interaction function in equation (2.3). If 2r 0 < d−ǫ, the quadratic term will be constant, and no force will exist to stabilize the circle. In order to use equation (3.6) then, we have to ensure that we are on the upper branch of figure 2.
Equation (3.6) gives the value of β C that provides an extremum of e 0 with respect to changes of radius a 0 at a givenr 0 (e 1 (r 0 ) = 0), but we still need to check that the circle of radiusr 0 is indeed stable to perturbations with non-zero frequency, i.e. that e 2 (k,r 0 ) is non-negative for all k. Scaling arguments mean that in fact the sign of e 2 depends only on the combinationsr 0 = r 0 /d andα C = (d/λ C )α C . The equation for e 2 can then be used to obtain bounds onα C in terms ofr 0 . (Details of these calculations and bounds can be found in [17].) The right hand side of figure 1 shows a plot of e 2 (k,r 0 ) againstr 0 k for the same parameter values used for the left hand side, showing that it is non-negative for allr 0 k.
We call the resulting model, the energy E g with parameters chosen according to the above criteria, the 'gas of circles' model.

Geometric experiments
To illustrate the behaviour of 'gas of circles' model, in this section we show the results of some experiments using E g (there are no image terms). Figure 3 shows the result of gradient descent using E g starting from various different initial regions.
(For details of the implementation of gradient descent for higher-order active contour energies using level set methods, see [14].) In the first column, four different initial regions are shown. The other three columns show the final regions, at conver-gence, for three different sets of parameters. In particular, the three columns havê r 0 = 15.0, 10.0, and 5.0 respectively.
In the first row, the initial shape is a circle of radius 32 pixels. The stable states, which can be seen in the other three columns, are circles with the desired radii in every case. In the second row, the initial region is composed of four circles of different radii. Depending on the value ofr 0 , some of these circles shrink and disappear. This behaviour can be explained by looking at figure 1. As already noted, the energy of a circle e 0 has a maximum at some radius r max . If an initial circle has a radius less than r max , it will 'slide down the energy slope' towards r 0 = 0, and disappear. If its radius is larger than r max , it will finish in the minimum, with radiuŝ r 0 . This is precisely what is observed in this second experiment. In the third row, the initial condition is composed of four squares. The squares evolve to circles of the appropriate radii. The fourth row has an initial condition composed of four differing shapes. The nature of the stable states depends on the relation between the stable radius,r 0 , and the size of the initial shapes. Ifr 0 is much smaller than an initial shape, this shape will 'decay' into several circles of radiusr 0 .

Likelihood energy and experiments
In this section, we apply the 'gas of circles' model to the extraction of trees from aerial images. We give a brief state of the art for tree crown extraction, and then present the likelihood energy we use in section 4.2. In section 4.3, we describe tree crown extraction experiments on aerial images and compare the results to those found using a classical active contour model. In section 4.4, we examine the robustness of the method to noise using synthetic images. This illuminates the principal failure modes of the model, which will be further discussed in section 5, and (Initial) compare the results to those obtained using a classical active contour model.

Previous work
The problem of locating, counting, or delineating individual trees in high resolution aerial images has been studied in a number of papers. For example, [18] observes that trees are brighter than the areas separating them. Local minima of the image are found using a 3 × 3 filter, and the 'valleys' connecting them are then found using a 5 × 5 filter. The tree crowns are subsequently delineated using a five-level rule-based method designed to find circular shapes, but with some small variations permitted. While the method is quite effective in separating trees, the size of the filters results in an significant overestimation of the size of the trees. [19] concentrates on spruce tree detection using a template matching method. The 3D shape of the tree is modelled using a generalized ellipsoid, while illumination is modelled The above methods use a series of ad hoc steps rather than a single unified model, which makes identifying the assumptions behind the method difficult. Closer in spirit to the present work is that of [21], which models the collection of tree crowns by a marked point process, where the marks are circles or ellipses. An energy is defined that penalizes, for example, overlapping shapes, and controls the parameters of the individual shapes. Compared to the work described in this paper, the method has the advantage that overlapping trees can be represented as two separate objects, but the disadvantage that the tree crowns are not precisely delineated due to the small number of degrees of freedom for each mark.

Likelihood energy and gradient descent
In order to couple the region model E g to image data, we need a likelihood, P(I|R, K).
The images we use for the experiments are coloured infrared (CIR) images. Originally they are composed of three bands, corresponding roughly to green, red, and near infrared (NIR). Analysis of the one-point statistics of the image in the region corresponding to trees and the image in the background, shows that the 'colour' information does not add a great deal of discriminating power compared to a 'greyscale' combination of the three bands, or indeed the NIR band on its own.
We therefore model the latter.
The image resolution is ∼ 0.5m/pixel, and tree crowns have diameters of the order of ten pixels. Little dependence remains between the pixels at this resolution, which means, when combined with the paucity of statistics within each tree crown, that pixel dependencies (i.e. texture) are very hard to use for modelling purposes. We therefore model the interior of tree crowns using a Gaussian distribution with mean µ and covariance σ 2 δ R , where δ A is the identity operator on images on A ⊂ Ω.
The background is very varied, and thus hard to model in a precise way. We use a Gaussian distribution with meanμ and varianceσ 2 δR. In general, µ >μ, and σ <σ; trees are brighter and more constant in intensity than the background. The boundary of each tree crown has significant inward-pointing image gradient, and although the Gaussian models should in principle take care of this, we have found in practice that it is useful to add a gradient term to the likelihood energy. Our likelihood thus has three factors: where I R and IR are the images restricted to R andR respectively, and g R and gR are proportional to the Gaussian distributions already described, i.e.
and similarly for gR. The function f ∂R depends on the gradient of the image ∂I on the boundary ∂R: where n is the unnormalized outward normal to γ. The normalization constant Z is thus a function of µ, σ,μ,σ, and λ i . Z is also a functional of the region R. To a first approximation, it is a linear combination of L(∂R) and A(R). It thus has the effect of changing the parameters λ C and α C in E g . However, since these parameters are essentially fixed by hand (the criteria described in section 3.1.1 only allow us to fix β C and constrain α C ), knowledge of the normalization constant does not change their values, and we ignore it once the likelihood parameters have been learnt.

The full model is then given by
The energy is minimized by gradient descent. The functional derivatives of all terms except the quadratic term in E g are standard. The functional derivative of the quadratic term gives rise to a gradient descent force given bŷ To evolve the region we use the level set framework of [22] extended to the demands of nonlocal forces such as equation (4.3) [14].
The computational complexity of the algorithm is unknown without a bound on the number of iterations. However, the complexity of one iteration is easily analysed. Equation

Tree crown extraction from aerial images
In this section, we present the results of the application of the above model to 50 In the experiments, we compare our model to a classical active contour model (β C = 0). The parameters µ, σ,μ, andσ were the same for both models, and were learned from hand-labelled examples in advance. The classical active contour prior model thus has three free parameters (λ i , λ C and α C ), while the 'gas of circles' model has six (λ i , λ C , α C , β C , d and r 0 ). We fixed r 0 based on our prior knowledge of tree crown size in the images, and d was then set equal to r 0 . Once α C and λ C have been fixed, β C is determined by equation (3.6). There are thus three effective parameters for the HOAC model. In the absence of any method to learn λ i , α C and λ C , they were fixed by hand to give the best results, as with most applications of active contour models. The values of λ i , λ C and α C were not the same for the classical active contour and HOAC models; they were chosen to give the best pos-          Table 1 shows the percentages of correct tree detections, false positives and false negatives (two joined trees count as one false negative), obtained with the classi-  Table 1 Results on real images using a classical active contour model (left) and the 'gas of circles' The typical runtime of the 'gas of circles' model in these experiments (image size O(100) pixels) is of the order of ten minutes on a normal personal computer. In our implementation, this is approximately ten times slower than using classical active contours.

Noisy synthetic images
In this section, we present the results of tests of the sensitivity of the model to noise in the image. Fifty synthetic images were created, each with ten circles with radius 8 pixels and ten circles with radius 3.5 pixels, placed at random but with overlaps rejected. Six different levels of white Gaussian noise, with image variance to noise power ratios from −5 dB to 20 dB, were then added to the images to generate 300 noisy images. Six of these, corresponding to noisy versions of the same original image, were used to learn µ, σ,μ, andσ. The model used was the same as that used for the aerial images, except that λ i was set equal to zero. The parameters were adjusted to give a stable radius of 8 pixels.
The results obtained on the noisy versions of one of the fifty images are shown in figure 8. Table 2 shows the proportion of false negative and false positive circle detections with respect to the total number of potentially correctly detectable circles (500 = 50 × 10), as well as the proportion of 'joined circles', when two circles are grouped together (an example can be seen in the bottom right image of figure 8).
Detections of one of the smaller circles (which only occurred a few times even at the highest noise level) were counted as false positives. The method is very robust with respect to all but the highest levels of noise. The first errors occur at 5 dB, where there is a 2% false positive rate. At 0 dB, the error rate is ∼ 10%, i.e. one of the ten circles in each image was misidentified on average. At −5 dB, the total error rate increases to ∼ 30%, rendering the method not very useful.
Note that the principal error modes of the model are false positives and joined circles. There are good reasons why these two types of error dominate. We will discuss them further in section 5.

Circle separation: comparison to classical active contours
In a final experiment, we simulated one of the most important causes of error in tree crown extraction, and examined the response of classical active contour and HOAC models to this situation. The errors, which involve joined circles similar to those found in the previous experiment, are caused by the fact that in many cases nearby tree crowns in an image are connected by regions of significant intensity with significant gradient with respect to the background, thus forming a dumbbell shape. Calling the bulbous extremities, the 'bells', and the join between them, the 'bar', the situation arises when the bells are brighter than the bar, while the bar is in turn brighter than the background, and most importantly, the gradient between the background and the bar is greater than that between the bar and the bells.
The first row of figure 9 shows a sequence of bells connected by bars. The intensity of the bar varies along the sequence, resulting in different gradient values. We applied the classical active contour and 'gas of circles' models to these images.  Table 2 Results on synthetic noisy images. FP, FN, J: percentages of false positive, false negative, and joined circle detections respectively, with respect to the potential total number of correct detections.
The middle row of figure 9 shows the best results obtained using the classical active contour model. The model was either unable to separate the individual circles, or the region completely vanished. The intuition is that if there is insufficient gradient to stop the region at the sides of the bar, then there will also be insufficient gradient to stop the region at the boundary between the bar and the bells, so that the region will vanish. On the other hand, if there is sufficient gradient between the bar and the background to stop the region, the circles will not be separated, and a 'bridge' will remain between the two circles. 2 The corresponding results using the 'gas of circles' model are shown in the bottom row of figure 9. All the circles were segmented correctly, independent of the gray level of the bar. Encouraging as this is, it is not the whole story, as we indicated in section 4.4. We make a further comment on this issue in section 5. Either the circles are not separated or the region vanishes. Bottom: the results using the 'gas of circles' model (2, 1, 5, 4.0, 8, 8). All the circles are segmented correctly.

Conclusion
Higher-order active contours allow the inclusion of sophisticated prior information in active contour models. HOACs are particularly well adapted to cases in which the topology is unknown a priori. In this paper, we have shown via a stability analysis that a HOAC energy can be constructed that describes a 'gas of circles', that is, it favours regions composed of an a priori unknown number of circles of a certain radius. The requirement that circles be stable, i.e. local minima of the energy, fixes one of the prior parameters and constrains the others.
The 'gas of circles' model has many uses in computer vision and image processing.
Combined with a suitable likelihood, we have applied it to the problem of tree crown extraction from aerial images of plantations. It performs better than simpler techniques such as maximum likelihood and classical active contours. In particular, it is better able to separate trees that appear joined in the data than is a classical active contour model.
The model is not without its issues, however. First, the computation time is too long. We are currently working on a phase field HOAC [23] version of the 'gas of circles' model that we hope will significantly reduce this time. Second, there are two significant error modes, as shown in the noise experiments of section 4.4: circles are found where the data does not ostensibly support them ('phantom circles'), and two circles may be joined into a dumbbell shape and never separated.
We discuss these in turn.
The first issue is that of 'phantom' circles. Circles of radiusr 0 are local minima of the prior energy. It is the effect of the data that converts such configurations into global minima. Were we able to find the global minimum of the energy, this would be fine. However, gradient descent finds only a local minimum. This can create problems in areas where the data does not support the existence of circles because a circle, once formed during gradient descent, cannot disappear unless there is an image force acting on it. We thus find that circles can appear and remain even though there is no data to support them.
The second issue is that of joined circles, discussed in section 4.5. Although the current HOAC model is better able to separate circles than a classical active contour, it still fails to do so in a number of cases, leaving a bridge between the circles. The issue here is a delicate balance between the parameters, which must be adjusted so that the sides of the bridge attract one another, thus breaking the bridge, and so that nearby circles repel one another at close range, so that the bridge does not re-form.
Again, this is at least in part an algorithmic issue. Even if the two separated circles have a lower energy than the joined circles, separation may never be achieved due to a local minimum caused by the bridge.
We propose to solve the first problem via a more detailed theoretical analysis of the circle energy that will allow us to remove the local minima causing the problem, and the second via an in-depth analysis of the energy of the dumbbell configuration.
Both these studies should lead to further constraints on the parameters, which is a desirable goal in itself.

Acknowledgements
This work was partially supported by EU project IMAVIS (FP5 IHP-

A.1 Linear terms
To compute the length, we need the magnitude of τ to second order. The metric in polar coordinates is ds 2 = dr 2 + r 2 dθ 2 , so we have that |τ (t)| 2 =ṙ(t) 2 + r(t) 2 by

equations (A.2). Substituting from equation (A.1) and (A.2) gives
Taking the square root, expanding it as √ 1 + x ≈ 1 + 1 2 x − 1 8 x 2 , and keeping terms to second order in the a k then gives Using equation (A.5), the boundary length is then given to second order by where we have used the reality of δr to set a −k = a * k , where * indicates complex conjugation, and orthonormality of the Fourier basis elements.
We can write the interior area of the region as Thus, using equations (A.1), and again using orthonormality, we have that

A.2 Quadratic terms
To compute the expansion of the quadratic term in equation (2.2) for E g , we need the expansions of τ (t) · τ (t ′ ) and Φ(R(t, t ′ )).
where unprimed quantities are evaluated at t and primed quantities at t ′ .

A.3 Combining terms
Using the above equations gives G(t, t ′ ) = r 2 0 cos(∆t)Φ(X 0 ) F 00 , even + (δr + δr ′ ) r 0 cos(∆t) Φ(X 0 ) + r 0 sin ∆t 2 Φ ′ (X 0 ) where the F ij denote the functions appearing in the terms of G, and 'odd' and 'even' refer to parity under exchange of t and t ′ . Each line, and hence G, is symmetric in t and t ′ , as it should be. We can now substitute the expressions for δr andδr in terms of their Fourier coefficients, and calculate π −π dt dt ′ G(t, t ′ ). We note that in the terms involving F 10 , F 11 , F 20 , F 22 , and F 23 , the presence of the symmetric or antisymmetric factors in δr and δr ′ simply leads to a doubling of the value of the integral for one of the terms in these factors, due to the corresponding symmetry or antisymmetry of the F functions. We therefore only need to evaluate one of these integrals for the relevant terms.