Amoeba Techniques for Shape and Texture Analysis
Abstract
Morphological amoebas are imageadaptive structuring elements for morphological and other local image filters introduced by Lerallut et al. Their construction is based on combining spatial distance with contrast information into an imagedependent metric. Amoeba filters show interesting parallels to image filtering methods based on partial differential equations (PDEs), which can be confirmed by asymptotic equivalence results. In computing amoebas, graph structures are generated that hold information about local image texture. This paper reviews and summarises the work of the author and his coauthors on morphological amoebas, particularly their relations to PDE filters and texture analysis. It presents some extensions and points out directions for future investigation on the subject.
Keywords: Adaptive morphology Morphological amoebas Curvaturebased PDE Selfsnakes Geodesic active contours Graph index Texture
1 Introduction
Mathematical morphology [38, 45, 46] has developed since the 1960s as a powerful theoretical framework from which versatile instruments for shape analysis in images can be derived, such as for structurepreserving denoising or shape simplification [23]. The fundamental building blocks of classical mathematical morphology are nonlinear local image filters like dilation, erosion, and median filters. They rely on aggregating intensities within a neighbourhood of any given pixel by e.g. maximum, minimum, and median operations. The selection of neighbourhoods for processing is classically done by shifting a sliding window of fixed size and shape across the image. In the context of morphology, this sliding window is known as structuring element.
More recently, concepts for adaptivity have been developed generally in image filtering and also specifically in morphology [51, 5]. One recent concept for adaptive morphology are morphological amoebas introduced by Lerallut et al. [32]. These are spacevariant structuring elements constructed from a combination of spatial distance measurement with local contrast measurement via an amoeba metric.
In earlier work by the author of the present paper and his coauthors, properties of amoeba filters and their relations to image filters based on partial differential equations (PDEs) were investigated [54, 57, 58]. As an application to image segmentation, an amoebabased active contour method was designed [53, 54, 56]. Recently, a combination of edgeweighted graphs generated in the computation of amoebas with graph indices was used to introduce a new class of texture descriptors [55] which are currently under further investigation. This paper reviews and summarises the results from these works. Directions of ongoing research on this topic are sketched.
With focus on giving a comprehensive overview of the theory that has been developed in various earlier publications, the (mostly lengthy) proofs of the results are omitted here and referred to the respective original sources. Nevertheless, the main principles underlying the proofs are shortly outlined. Although amoeba filtering of multichannel images has been addressed to some extent in [57], this aspect of the topic presents itself in a stage too early for a summarised presentation, and is therefore not included in the present paper.
In the following the structure of the paper is detailed, highlighting contributions that are novel in this presentation.
Section 2 introduces the concept of morphological amoebas as imageadaptive structuring elements in the spacediscrete as well as the spacecontinuous setting. To ease bridging to the graph techniques discussed later in Section 6, the presentation in the discrete case emphasises the modelling of discrete images by neighbourhood graphs and uses standard terminology from graph theory, thereby following [55]. The presentation of the spacecontinuous case is similar to that e.g. in [57].
The application of amoebas in image filtering is the topic of Section 3. Median filters, morphological dilation and erosion are presented together with their relationship to PDE image filters, reproducing herein results from [54, 56, 57, 58]. Regarding the association between amoeba metrics on the discrete filtering side and edgestopping functions occurring in the corresponding PDEs, the current work adds to the previously considered exemplary and (Euclidean) amoeba metrics as a third simple case the (maximum) amoeba metric and states explicitly the corresponding edgestopping function. Moreover, the amoeba variants of morphological opening and closing are included in the description for the first time. For dilation, erosion, opening and closing filters, the presentation here emphasises the algebraic background including maxplus/minplus convolution and conjugacy of structure elements.
Section 4 considers the application of amoeba techniques to devise basic algorithms for unsupervised segmentation of greyvalue images, namely the amoeba active contours (AAC) first introduced in [53] and further investigated in [54, 56]. Results from [56] on the relation between AAC and geodesic active contours are reported.
In image filtering by nonlinear PDEs, one often computes the nonlinearities not from the input images themselves but from Gaussian presmoothed versions of these, in order to reduce noise sensitivity of filters and to improve numerical stability. This is also the case with selfsnakes and active contour PDEs; note that the selfsnakes PDE is even illposed without such presmoothing. Section 5 investigates the effect of presmoothing in the selfsnakes PDE using perturbation analysis on a synthetic example; furthermore, it discusses how a comparable stabilisation can be achieved in the amoeba median filter framework. The analysis presented in this section relies on previous work in [54, 57] in which oscillatory perturbations aligned with the gradient direction were studied, and extends it by including also perturbations aligned with the level line direction.
Section 6 is devoted to a different direction of application of amoeba ideas. Noticing that the computation of discrete amoeba structuring elements is intimately related with graph structures – a weighted neighbourhood graph, weighted and unweighted Dijkstra search trees – in the neighbourhood of each pixel, one can try to extract local texture information from these graphs. Quantitative graph theory [13] offers a variety of graph indices for generating quantitative information from graph structures. The presentation of the construction of texture descriptors from amoebas and graph indices in this section follows [55]. Compared to the large set of descriptors covered in [55], only a few representatives are shown here, complementing their mathematical description by a visualised example. Extending the previous work on texture discrimination in [55], the present paper also shows a first example of the new texture descriptors in texture segmentation by using the descriptors as components of an input image for multichannel GAC segmentation.
2 Morphological Amoebas
Wellknown local image filters such as the mean filter, median filter, morphological dilation or erosion consist of two steps: a slidingwindow selection step, and the aggregation of selected input values by taking e.g. the arithmetic mean, median, maximum or minimum. A strategy to improve the sensitivity of such filters to important image structures is to modify the selection step by using spatially adaptive neighbourhoods instead of a fixed sliding window. The general idea is to give preference in the selection to neighbouring image locations with similar intensities, and thus to reduce the flow of greyvalue information across high contrast steps or slopes in the filter process.
First introduced by Lerallut et al. [32, 33] as structuring elements for adaptive morphology, morphological amoebas are a specific type of such spatially adaptive neighbourhoods. Their construction relies on the combination of spatial distance in the image domain with greyvalue contrast into a modified metric on the image.
2.1 Edgeweighted Neighbourhood Graph
To define morphological amoebas on discrete images, we start by considering edgeweighted graphs based on the image grid.
Definition 1.
Let be a discrete image. Construct an edgeweighted graph with vertex set , edge set and weights as follows. The vertex set is formed by all pixels of . Two vertices , are connected, , if and only if pixels , are neighbours under a suitably chosen neighbourhood notion. To define the edge weights for an edge , consider the corresponding pixel locations and as well as the intensities and , and set to
(1) 
where denotes Euclidean distance in the image plane, is a contrast scale parameter weighting between spatial and tonal distances, and is a norm on which can be rewritten as
(2) 
with a monotonically increasing function (by continuity, ).
The edgeweighted graph is called neighbourhood graph of .
In this definition, neighbourhood can be understood as a 4neighbourhood, as done in [32], or as an 8neighbourhood as in [55, 57, 58]. The latter choice gets somewhat closer to a Euclidean measurement of spatial distances in the image plane and is therefore also considered the default in the present work.
As to the norm function , the setting corresponds to the metric also used in [32] that gives
(3) 
whereas entails a Euclidean () metric in which the edge weights are obtained by the Pythagorean sum
(4) 
A straightforward generalisation is
(5) 
which in the limit also includes and the corresponding edge weight
(6) 
2.2 Discrete Amoeba Metric
We use the edgeweighted neighbourhood graph to define the discrete amoeba metric on image .
Definition 2.
Let a discrete image be given. Let be its neighbourhood graph with edge weights given by (1). Define for two pixels and their distance as the minimal total weight (length) among all paths between and in . Then is called (discrete) amoeba metric on .
The metric is called amoeba metric, , if it is derived from (5), or amoeba metric if it is obtained from . The amoeba metric is also called Euclidean amoeba metric.
Definition 3.
In a discrete image with amoeba metric , an amoeba structuring element (short: amoeba) with amoeba radius and reference point at pixel is a discrete ball around pixel in the amoeba metric, i.e. the set of all vertices within a distance from ,
(7) 
The derivation of amoebas from a metric with a global radius parameter has an interesting consequence: for two pixels , , one has
(8) 
which is helpful in the design of some morphological filters.
2.3 Computation of Discrete Amoebas
To compute amoebas in a discrete image, one has to search the neighbourhood of each given reference pixel in order to identify the pixels with amoeba distance . Given that the edge weights in are nonnegative, this can be achieved by running Dijkstra’s shortest path algorithm [16] on starting at pixel . As this algorithm enumerates neighbour pixels in order of increasing path weight, it can be stopped as soon as a pixel with is visited.
Moreover, by the construction of the amoeba distance it is clear that the Euclidean distance in the image domain is a lower bound for the amoeba distance between pixels. Therefore the Dijkstra algorithm for the start vertex can be run on the subgraph of that contains just the pixels from the Euclidean neighbourhood of .
2.4 Amoebas on Continuous Domains
Even superficial inspection of results obtained by some amoeba filters indicates that they have striking similarities to image processing methods based on partial differential equations (PDEs). This observation has been substantiated in [56, 57, 58] by studying spacecontinuous versions of amoeba filters; the results proven there allow to interpret amoeba filters as time steps of explicit discretisations for suitable PDEs.
To devise spacecontinuous versions of amoeba filters, one has to translate first the notion of amoeba metric to the spacecontinuous setting. Once this is done, the definition of an amoeba as a ball around a reference point is straightforward.
The amoeba metric for a spacecontinuous greyvalue image – a realvalued function over a connected compact image domain – can be stated by assigning to each two given points as their distance the minimum of a path integral between and . Just like the edge weights in the discrete amoeba construction, the integrand of the path integral is obtained by applying a suitable norm to the spatial metric (the Euclidean curve element of the path) and the greyvalue metric (the standard metric on the real domain), such that the amoeba distance reads as
(9) 
where runs over all regular curves with , , and can be chosen as in the discrete case.
Let us associate to the function its (vertically rescaled) graph, the manifold . Then we see that the amoeba distance between two points , in the image domain can be interpreted as a distance on . The points herein are given by , . To define the metric on , consider a metric in the surrounding space that combines Euclidean metric in the plane with the standard metric in direction via the function from (2) that appeared already in the original construction of the amoeba metric. Using in , the metric is obtained as its induced metric on the submanifold . Figure 1 illustrates that the amoeba structuring element is then the projection of a unit disk on back to the image plane.
3 AmoebaBased Image Filters
To obtain applicable image filters, the amoeba procedure described above is used as a selection step and needs to be complemented by some aggregation step. We consider here standard choices of aggregation operators from classical local filters; introducing also modifications into this part of the filtering procedure is left as a possible direction for future research. Moreover, keeping close to the original context in which amoebas were developed, we focus on morphological operators. Here, morphological operators are characterised by their invariance under arbitrary monotonically increasing transformations of the intensities, see e.g. [37], which means that also median and quantiles belong to this class.
3.1 Median
A median filter aggregates the intensity values of the selected pixels by taking their median. In the nonadaptive, slidingwindow setting this filter can be traced back to Tukey [50], and since then it has gained high popularity as a simple denoising filter that preserves discontinuities (edges) and its robustness with respect to some types of noise. Median filtering can be iterated. Unlike average filters, the median filter on a discrete image possesses nontrivial steady states, socalled root signals [17], that depend on the filter window. The smaller the filter window, the faster the iterated median filtering process locks in at a root signal.
Despite the nice preservation of edges, the nonadaptive median filter involves a displacement of curved edges in inward direction and rounding of corners that is often undesired. Amoeba median filtering greatly reduces this effect. Figure 3 demonstrates this by an example.
PDE Approximation
As noticed already in 1997 by Guichard and Morel [20], the overall robust denoising effect and the characteristic cornerrounding behaviour of standard median filtering resemble the properties of the wellknown (mean) curvature motion PDE [1]. Further analysis confirmed this observation by proving an asymptotic relationship between the two filters, as set forth in the following proposition.
Proposition 1 (Guichard and Morel, 1997 [20]).
For a smooth function , one iteration of median filtering with a ball as structuring element approximates for a time step of size of the curvature motion PDE [1]
(10) 
This seminal result motivates the investigation of relations between amoeba and PDE filters whose results are reviewed in the further course of the present paper.
Just like amoeba median filtering differs from standard median filtering by an adaptation procedure that suppresses smoothing across edges, the curvature motion equation (10) has a counterpart in which also the flow across edges is suppressed. This socalled selfsnakes filter [44] allows curvaturebased image smoothing and simplification, preserves and even enhances edges, while at the same time avoiding to shift them, as curvature motion does. It turns out that indeed amoeba median filtering is connected to selfsnakes by a similar asymptotic relationship as that of Proposition 1, as follows.
Theorem 1 ([57, 58]).
For a smooth function , one iteration of amoeba median filtering with amoeba radius approximates for a time step of size of the selfsnakes PDE [44]
(11) 
where is a decreasing edgestopping function that depends on the amoeba metric being used.
Proofs for Theorem 1 have been given in [57, 58]. While these proofs are not reproduced in detail here, it is of interest to describe the two different strategies that are used in these proofs. These approaches form also the basis for the further amoeba—PDE asymptotics results presented in Section 3.2.
Proof Strategies
The crucial observation for all median filter—PDE equivalence results since Guichard and Morel’s proof of Proposition 1 in [20] is that the median of a smooth function within a given compact structuring element is the function value whose corresponding level line divides the structuring element into two parts of equal area. Herein it is assumed that each value of within the structuring element is associated with a unique level line segment inside , which is satisfied for sufficiently small fixed or amoeba structuring elements whose reference point is not an extremum of , and therefore acceptable when studying the limit .
The amount by which a single median filtering step changes the function value at the reference point of the structuring element then corresponds, up to multiplication with , to the distance between the areabisecting level line and the level line through , see the illustration in Figure 4(a). The two approaches discussed in the following differ in the way how they measure the area of the structuring elements and parts thereof.
Proof Strategy I. The first strategy has been followed in [58] to prove Theorem 1 for the entire class of amoeba metrics discussed in Section 2 above, see also the more detailed version in [57, Section 4.1.1]. It is close to the approach from [20] in that it develops the smooth function around the reference point into a Taylor expansion up to second order. The Taylor expansion is then used to approximate, for an amoeba , three items: first, the range of function values occurring within , i.e. the minimum and maximum , second, the length of the level line segment for each , and third, the density of level lines around each , which equals the steepness of the slope of near the level line of .
Integrating the lengths of level lines over function values, weighted with their reciprocal densities, yields the area of , i.e.
(12) 
As this integral effectively runs over level lines, splitting the integration interval exactly corresponds to cutting at some level line. The calculation of the desired median of within is then achieved by determining a suitable splitting point in the integration interval so that the integrals on both subintervals become equal.
Summarising, this strategy describes the amoeba shape in terms of a curvilinear coordinate system aligned with the gradient and level line directions at , in which the level lines take the role of coordinate lines, compare Figure 4(b).
Proof Strategy II. The second strategy abandons the consideration of the individual level lines within ; the only level line that is explicitly studied is the one through itself. Instead of the distorted Cartesian coordinate system one uses polar coordinates to describe the shape of the amoeba. This approach has first been used in [54] in the context of amoeba active contours (see Section 4.1), and again in [57, Section 4.1.2], both times restricted to the Euclidean amoeba metric. It has been extended to cover the full generality of amoeba metrics under consideration in [56], again for amoeba active contours.
Writing the outline of as a function of the polar angle , the amoeba’s area is stated by the standard integral for areas enclosed by function graphs in polar coordinates as
(13) 
Unlike for (12), splitting this integral yields areas of sectors instead of segments; however, if the level line through happens to be a straight line, splitting up the integral (13) at the pair of opposite angles corresponding to the level line direction yields the areas of two segments into which is cut by that level line, compare Figure 5(a).
Provided that is symmetric (w.r.t. point reflection at the reference point), the two segments are of equal area, making in this case the median equal to . Deviations from this situation that make the median differ from can be separated into two contributions: first, the asymmetry of the amoeba; second, the curvature of the level lines. Crosseffects of the two contributions influence only higher order terms that can be neglected in the asymptotic analysis; thus the two sources can be studied independently. In approximating the area difference caused by the asymmetric amoeba shape, one can assume that the level lines are straight, see Figure 5(b), while the level line curvature effect can be studied under the assumption that has symmetric shape, see Figure 5(c).
Finally, the combined effect must be compensated by a parallel shift of the level line through , compare again Figure 4(a). From the shift the median, and thus the righthand side of the PDE approximated by the amoeba filter, can be derived.
Amoeba Metrics and EdgeStopping Functions
It remains to specify the relation between amoeba metric and edgestopping function mentioned in Theorem 1. In [58, 57], the following representation of in terms of the function defining the amoeba metric has been proven.
(14) 
where is short for , i.e. the square of the inverse of , and for the cube .
In the case of the Euclidean amoeba metric, , the expression (14) simplifies to
(15) 
which is, up to the substitution , the PeronaMalik diffusivity [39] that is also one of the common choices for in the selfsnakes equation.
When using the amoeba metric, , the integral in (14) can be numerically evaluated, and one obtains an edgestopping function that differs from (15) in that it decreases away from already with nonvanishing negative slope, thus reacting more sensitive to even small image contrasts.
Finally, for the amoeba metric, , it is again possible to state in closed form,
(16) 
which shows that is completely insensitive to image contrasts up to and then starts decreasing with a kink. All three edgestopping functions are depicted in Figure 6.
3.2 Dilation and Erosion
The two most fundamental operations of mathematical morphology, dilation and erosion, use as aggregation step the maximum and minimum of intensities, respectively. This can naturally be done also in combination with an amoebabased pixel selection step.
We point out that the standard dilation of an image with fixed structuring element can be written as
(17) 
where denotes the function
(18) 
The last term in (17) allows an interesting interpretation in terms of the maxplus algebra [3, 42], an algebraic structure on in which the maximum operation takes the role of addition in the usual algebra of real numbers, while addition takes the role of multiplication. It is evident that (17) is nothing else but a convolution of and in the maxplus algebra, see [36].
In writing erosion in an analogous way, we follow a convention frequently used in the literature by using instead of the structuring element the conjugate structuring element , which comes down geometrically to a point reflection on the origin, . The advantage of this convention is that subsequent definitions like those for opening and closing become simpler [24], compare Section 3.3.
Defining then as zero on , but outside, erosion is stated as
(19) 
which can be interpreted again as a convolution of and in the minplus algebra [36].
Abandoning the fixed window and using a family of structuring elements located at pixel , one can write amoeba dilation as
(20)  
(21) 
Just as the last term in (17) is a maxplus convolution, the righthand side (20) is the maxplus analogon of a (discretised) integral operator. Herein, acts as the maxplus counterpart of just the same type of integral kernel that appears as pointspread function in spacevariant image deconvolution models.
Similarly, amoeba erosion becomes a minplus integral operator with a minplus kernel . Generally, conjugate structuring elements in the spacevariant case are given by
(22) 
Interestingly, if is made up by amoebas , there is no difference whether the conjugate structuring elements or standard structuring elements are used in erosion: property (8) of the amoebas entails for all , or equivalently
(23) 
We will denote this property as selfconjugacy of amoebas.
Figure 8 shows the results of nonadaptive and amoeba dilation and erosion of an example image depicted in Figure 7. Nonadaptive as well as amoebabased dilation extend bright image details, but it can be seen that the spreading of bright image parts is stopped at strong edges; similarly for the propagation of dark details by erosion.
PDE Approximation
It is a wellknown fact that HamiltonJacobi PDEs
(24) 
describe dilation (“” case) and erosion (“”) of continuousscale images or levelset functions in the sense that evolution of an initial image by (24) up to time yields the dilation or erosion of with a Euclidean ballshaped structuring element of radius . It can therefore be expected that amoeba dilation and erosion, too, should be related to hyperbolic PDEs resembling (24). The following result from [57] confirms this intuition.
Theorem 2 ([57]).
For a smooth function , one step of amoeba dilation or amoeba erosion with amoeba radius and Euclidean amoeba metric approximates for a time step of size of an explicit time discretisation of the HamiltonJacobitype PDE
(25) 
where the “” sign applies for dilation, and “” for erosion.
Note that unlike in Theorem 1 the time step size here depends linearly, not quadratically, on . In [57] the theorem is formulated slightly more general to cover also amoeba quantile filters that interpolate in a natural way between median filtering (), dilation () and erosion (). As a result of the different order of decay of for , it comes as no surprise that for always the advection behaviour of the HamiltonJacobi equation (25) dominates over the parabolic equation (10), thus turning quantile filters into “slower” approximations to the same PDE.
3.3 Opening and Closing
In mathematical morphology, the opening of an image with (fixed) structuring element is defined as the concatenation of an erosion followed by a dilation with . In case is not pointsymmetric it is essential that, as mentioned in Section 3.2, the conjugate structuring element is used in the erosion step. Opening therefore reads as
(26) 
Analogously, closing is defined as dilation followed by erosion,
(27) 
Again, it is straightforward to turn these operations into adaptive variants by using amoeba structuring elements. Amoeba opening and closing of image with amoebas of radius are given as
(28)  
(29) 
where .
It is worth noticing that the difficulty about using the conjugate set of structuring elements for erosion disappears here due to the selfconjugacy (23) of the amoeba structuring element set.
As it is essential to use the same set of structuring elements in the dilation and erosion step, both steps must be carried out with the amoebas obtained from the original image. The underlying principle is that in the second step (dilation for opening or erosion for closing) each pixel should influence exactly those pixels which have influenced it in the first step before. As a consequence, e.g. amoeba opening is not exactly the same as amoeba erosion followed by amoeba dilation – this sequence would be understood by default as recalculating amoebas after the erosion step, i.e.
(30) 
which is inappropriate for an opening operation.
In Figure 9 exemplary results of nonadaptive and amoebabased closing and opening of the test image from Figure 7 are shown. Like its nonadaptive counterparts, amoebabased closing and opening remove smallscale dark or bright details, respectively. However, the amoeba versions do this in a less aggressive way. Extended narrow structures that are often removed partially by the nonadaptive filters are more often preserved as a whole, with reduced contrast, or removed completely by the amoeba filters, see e.g. the roof front edge descending to the right from the chimney, and the acute roof corner separating it from the sky.
Opening and Closing Scale Spaces and PDEs
The association between median, dilation and erosion filters and PDEs is inherently related to the scale space structures of these filters, compare [25]. All of these filters form an additive semigroup in the sense that iterative application of the same filter yields an increasing filter effect that naturally adds up over iteration numbers. In the case of dilation and erosion iteration numbers are also in linear relation with increasing structuring element size, as dilating an initial image times with (nonadaptive) structuring element radius is equivalent to dilating once with radius . Such an additive semigroup structure perfectly matches initial value problems for PDEs in which, too, evolution times add up.
While opening and closing, too, have a scale space structure, their semigroup operation is not additive but supremal as it is based on taking the maximum of parameters. For example, repeating the same opening or closing operation on a given image just reproduces the result of the first application of the filter (i.e., opening and closing operators are idempotent); and concatenating two openings or two closings with structuring element radii , gives an opening or closing with radius .
For this reason, also amoeba opening and closing are not associated with PDE evolutions in the same way as the previous filters. Possible relations to PDEbased filters may be considered in future research.
4 GreyScale Segmentation
Following established terminology, image segmentation denotes the task to decompose a given image into regions that are in the one or other way homogeneous in themselves but different from each other, with the intention that these regions are meaningful in that they are associated to objects being depicted. Intensitybased segmentation uses intensity as the main criterion of homogeneity within and dissimilarity between segments. Specialising to the case of two segments (foreground and background) with the additional geometric hypothesis that segments are separated by sharp and smooth contours, contourbased segmentation approaches based on curve or level set evolutions lend themselves as tools for segmentation, with active contours as an important representative. In this section we show how amoeba algorithms can be made useful in this context.
Despite the fact that experiments on magnetic resonance data are used to illustrate the concepts in this section, this is not meant to make a claim that neither active contour nor the related active region methods (which are not discussed further here) in their pure form could serve as a stateoftheart segmentation method for medical images. In fact, competitive results in medical image segmentation are nowadays achieved by complex frameworks that often include active contours and/or active regions as a component but in combination with additional techniques that allow to bring in anatomical knowledge such as shape and appearance models [11]. An early representative of these frameworks is [34], which has been followed by many more since then. Like classical geodesic active contours, the amoeba active contours presented in the following could be integrated into this type of framework but this has not been done so far.
4.1 Amoeba Active Contours
The standard procedure of an active contour, or snake, method starts with some initial contour which may be obtained automatically from some previous knowledge or heuristics regarding the position of a sought structure, or from human operator input. Representing this contour either by a sampled curve or by a levelset function, it is then evolved up to a given evolution time or up to a steady state by the action of some parabolic PDE, which is often derived as a gradient descent of a segmentation energy in the image plane. An important representative are geodesic active contours (GAC) [9, 30]. Their segmentation energy is essentially a curve length measure of the contour in a modified metric on the image plane that favours placing the contour in highcontrast locations. The PDE for GAC in levelset representation reads
(31) 
Herein, is the evolving levelset function in the plane that represents the actual evolving contour as one of its level sets (by default, the zerolevel set), and is the invariable image being segmented. The similarity of (31) to selfsnakes (11) (which were actually inspired from active contours, thus the name) together with the link between amoeba median filtering and self snakes established by Theorem 1 suggest that an amoeba median approach could be used to evolve the level set function instead of equation (31).
Introduced in [53], the resulting amoeba active contour (AAC) algorithm proceeds as follows:

Compute amoeba structuring elements based on the input image .

Initialise the evolving levelset function to represent the initial contour.

Evolve the image by median filtering with the amoebas from Step 1 as structuring elements.
Results from this algorithm look qualitatively fairly similar to those from GAC, as will also be demonstrated later in this section.
4.2 PDE Approximation
In order to study the relation between AAC and GAC, it makes sense again to consider a spacecontinuous model and to investigate the PDE approximated by AAC in the case of vanishing amoeba radius. The following result was proven in [56]. Note that in this theorem the contrast scale parameter is fixed to for simplicity, which, however, is no restriction of the result because in the active contour setting in question, the case is easily mapped to by just scaling the intensities of image by .
Theorem 3 ([56]).
Let a smooth levelset function be filtered by amoeba median filtering, where the amoebas are generated from a smooth image . Assume that the amoeba metric is given by (9), (2) with . One step of this filter for then approximates for a time step size of size of an explicit time discretisation of the PDE
(32) 
with the coefficients given by
(33)  
(34) 
Here, and are unit vectors in gradient and level line direction, respectively, for , whereas and are the corresponding unit vectors for , and is the angle between both gradient directions.
The proof of this result is found for the case of the Euclidean amoeba metric in [54], and for general amoeba metric in [56]. It relies on Proof Strategy II from Section 3.1.2.
An attempt to analyse AAC using Proof Strategy I had been made in [53], where, however, only a special case was successfully treated: The theorem proven in [53] states that AAC approximates the GAC equation (31) if image and level set function are rotationally symmetric about the same centre.
In fact, the rotational symmetry hypothesis can be weakened; what is needed for (32)–(34) to reduce to the exact GAC equation is actually, whenever (thus, , ), and hold, (32)–(34) boil down to the GAC equation (31).
At first glance, this is still a very artificial choice; however, looking at the geometrical implications of this setting, one sees that it means that the level lines of are aligned to those of , have the same curvature, and the image contrast in both and does not change along these level lines. Thereby the hypothesis of this special case is well approximated in the nearconvergence stage of a segmentation process when the object—background contrast is more or less uniform along the contour.
As a consequence, the coincidence of AAC and GAC in this case justifies that both approaches can expected to yield very similar types of segmentations. The convergence behaviour towards these segmentations may differ more; a closer comparison of both PDEs in [54, 56] based on typical geometric configurations indicates that the amoeba active contour PDE drives contours toward image contours in a more pronounced way.
Figure 10 presents an example that confirms the overall similarity between amoeba and geodesic active contours but also the tendency of AAC to adapt more precise to very smallscaled edge details. Frame (a) shows the original image with an initial contour roughly enclosing the cerebellum. Frames (b) and (c) demonstrate segmentation by AAC with Euclidean and amoeba metrics, respectively, while Frame (d) shows a GAC result for comparison.
4.3 Force Terms
Geodesic active contours in their basic form (31) suffer from some limitations. First of all, when initialised with a contour enclosing a large area with one or several small objects inside, the active contour process spends plenty of evolution time to slowly move the contour inwards until it hits an object boundary, due to the initially small curvature of the contour. Secondly, for pronounced concave object geometries, the process tends to lock in at undesired local minima that detect well some convex contour parts but shortcut concave parts via straight line segments. Similar problems can occur when segmenting multiple objects within one initial contour, see the examples in [31]. Thirdly, as the basic curvature motion process involves only inward movement of contours, it is generally not possible with (31) to segment objects from initial contours inside the object, which is sometimes desirable in applications. Due to their similarity to GAC, amoeba active contours share these problems.
A common remedy for these problems in the literature on active contour segmentation is the introduction of a force term. Its typical form is , i.e. essentially the righthand side of a HamiltonJacobi PDE for dilation or erosion, compare (24). An erosion force accelerates the inward motion of the contour; it allows to get past homogeneous areas faster, and helps the contour to find concave object boundaries and to separate multiple objects. By a dilation force it is possible to push the contour evolution in outward direction, which makes it possible to use initial contours inside objects.
In both cases, however, the force strength needs careful adjustment because dilation or erosion may also push the contour evolution across object boundaries, thereby preventing their detection.
In [10] where this modification was proposed first (by the name of “balloon force”), was chosen as constant, but the possibility to steer it contrastdependent, was mentioned. This has been done in [9, 31, 35] by modulating the force term in a geodesic active contour model with the same edgestopping function , such that the entire force term reads as with constant .
The relation between amoeba quantile filters and HamiltonJacobi PDEs mentioned in Section 3.2 indicates how to achieve a similar modification in the amoeba active contour algorithm: the median filter step should be biased, basically by replacing the median with some quantile. The most obvious way to do this is to use the quantile with a fixed . Within a discrete amoeba containing pixels, this means to choose the value ranked in the ordered sequence of intensities. However, taking into account that the amoeba size (or the amoeba area in the continuous setting) varies even for fixed with local image contrast, it is not less natural to think of as varying with the amoeba size. If one chooses inversely proportional to the amoeba size, this comes down to modify the median with a fixed rank offset , such that in an amoeba of pixels one would choose the intensity value with rank . These two variants of the AAC algorithm have been proposed in [53]. In [56] a third variant (“quadratic bias”) was introduced which chooses from the rank order the element with index with fixed . For these three scenarios, further analysis was provided in [56], based on the Euclidean amoeba metric. We summarise the results here.
Fixed offset bias. Choosing the entry at position from the rank order approximates a force term with . Note that in the symmetric case in which the PDE approximated by AAC coincides with the GAC equation this becomes exactly the “balloon force” term with constant dilation/erosion weight from [10].
Quadratic bias. Choosing the entry at index from the rank order of intensities yields an approximated force term . In the rotationally symmetric case this corresponds to the weight from [31].
To illustrate amoeba active contours with bias, Fig. 11 presents an example (shortened from [56]). Frame (a) is a test image with initial contour inside a mostly homogeneous object (the corpus callosum). Fig. 11(b) and (c) then show contours computed by amoeba active contours with fixed offset bias for two different evolution times, one intermediate, one displaying the final segmentation. For comparison, a segmentation with geodesic active contours is shown in (d).
We remark that in the AAC examples, a few pixels within the corpus callosum region are excluded from the segment, see the small isolated contour loops there. This is not a numerical artifact but a result from the precise adaption of amoebas to image structures even up to the resolution limit (pixel precision) of the image – the pixels not included in the segment are noise pixels with intensities significantly deviating from the neighbourhood, which are simply not included in any amoeba of outside pixels. Modifications like presmoothing input images can be applied to avoid this. On the contrary, the absence of such difficulties in the GAC example is a beneficial effect of the otherwise often undesirable numerical blurring effect of the finitedifference scheme.
5 PreSmoothing in SelfSnakes and Amoeba Filters
The approximation result of Theorem 1 associates iterated amoeba median filtering with the selfsnakes equation (11). Unlike (mean) curvature motion (10), selfsnakes possess edgeenhancing properties. Rewriting (11) by the product rule, one can state the selfsnakes process as
(35) 
in which the first summand is just a curvature motion process modulated by , whereas the second, advective, term is responsible for the edgeenhancing behaviour. Unfortunately, this term has a shockfilter property which makes not only its numerical treatment difficult – in finite difference schemes usually an upwind discretisation will be required to approximate it – but even entails illposedness of the PDE itself that is reflected in a noticeable staircasing behaviour. Indeed, as demonstrated by an experiment in [58], the result of a numerical computation of a selfsnakes evolution differs significantly if the underlying grid resolution is changed.
A common remedy to this illposed behaviour is to use presmoothing in the argument of the edgestopping function, i.e. to replace in (11) or (35) by where is the result of convolving with a Gaussian of standard deviation . Thereby, the illposedness of selfsnakes is removed, and a stable filtering achieved, at the cost of the additional smoothingscale parameter .
In this section, we deal with the question whether this staircasing phenomenon has also an analogue in the amoeba median filtering context, and what is an appropriate counterpart for the presmoothing modification on the amoeba side. This is done by quantitative analysis of a synthetic example, the first part of which has been published before in [54, 57].
5.1 PreSmoothing in Amoeba Median Filtering, and Amoeba Radius
First of all, notice that a straightforward translation of the presmoothing procedure to the amoeba median filtering context is to use in place of when computing the structuring elements in an amoeba median filtering step. This is actually an instance of the generalised amoeba median filtering procedure of the amoeba active contour setting, Sections 4.1 and 4.2, such that the PDE approximation result from Theorem 3 can be applied to see that it would approximate a PDE which is not identical to the standard selfsnakes with presmoothing, but closely related to it.
At second glance, however, it can be questioned whether the introduction of the smoothingscale parameter into the amoeba median filter is necessary. Unlike finitedifference schemes for selfsnakes, amoeba filtering by construction already involves a very similar smoothingscale parameter, namely, the amoeba radius . One can conjecture that the positive necessarily used in any amoeba computation could already provide a presmoothing effect similar to the Gaussian convolution in the PDE setting. This conjecture will be investigated in the following.
5.2 Perturbation Analysis of Test Cases
The starting point for constructing the test cases is a simple slope function that would be stationary under both selfsnakes and amoeba median filter evolutions, see Figure 12(a). From this slope, described by the function , , test cases are derived by adding small singlefrequency oscillations such as with frequency vectors .
Given the nonlinear nature of the filters under investigation, there is no superposition property for the effects of different perturbations of . Nevertheless, interactions between and the perturbations are always of higher order , such that the analysis of the firstorder effects of perturbations still gives a useful intuition about the behaviour of the filters.
Test Case 1: GradientAligned Oscillation
For the first test case, see [54, 57], the perturbation frequency is aligned with the gradient direction, , yielding the input signal schematically depicted in Figure 12(b),
(36) 
SelfSnakes Analysis. To determine the response of the selfsnakes evolution (35) to the perturbed signal (36), notice first that level lines of (36) are straight and parallel, such that one has and . Further, one has and , finally turning (35) into
(37) 
From this it can be read off that a frequency response factor occurs that grows indefinitely for high frequencies. Since the nonlinearity of (35) instantaneously spreads out the single perturbation frequency to higher harmonics, arbitrarily high amplification appears already within short evolution time, and the regularity of the evolving function is lost. This explains the staircasing behaviour of selfsnakes without presmoothing.
Using presmoothed in the edgestopping function argument, one has instead , and therefore
(38) 
with the frequency response factor that is globally bounded with its maximum at . Therefore, presmoothing ensures that the regularity of the evolving function is maintained.
Amoeba Filter Analysis. To analyse the effect of amoeba median filtering (with Euclidean amoeba metric) on the function (36), consider an amoeba of amoeba radius around , and assume that the contrast scale is chosen as .
The median of within that amoeba can be expressed via an integral formula, see [54, 57], which can be numerically evaluated to be approximately equal to with a frequency response factor . In other words, one amoeba median filter step amplifies the perturbation of (36) versus by the amplification factor .
Figure 13 shows results of numerical approximation of one amoeba median filtering step with , , on test images of type (36) with two different frequencies . The numerical computation was carried out on a discrete grid with mesh size . For best approximation to the spacecontinuous case, amoeba distances between pixels were computed by numerical integration instead of the Dijkstra search on the pixel graph. Denoting the filtered image by , numerical amplification factors can be computed as (with the usual scalar product of functions on a suitable bounded interval); these are in good accordance with the theoretical result.
Figure 14 shows the amplification function for together with its counterpart for one time step of selfsnakes with presmoothing, with the time step size matching the amoeba radius according to Theorem 1. The figure also includes numerical amplification factors for amoeba median filtering with the same parameters for frequencies . The parameter in the selfsnakes case has been chosen for a good match to the first wave of . With this parameter, the amplification behaviour for frequencies up to approx. is very similar for the presmoothed selfsnakes equation and amoeba median filtering. However, for higher frequencies the amplification factor of presmoothed selfsnakes rapidly approaches one (no amplification) whereas it oscillates around for the amoeba filter.
As a result, oscillations with sufficiently high frequency are just almost not amplified in the presmoothed selfsnakes evolution. With amoeba median filtering, they are amplified by the globally bounded factor in each iteration step. Whatever was in the initial image from (36), after a finite number of iterations the oscillations grow to a level for which the hypothesis of our analysis is no longer valid. Even in the spacecontinuous setting under consideration, oscillations cannot actually grow infinitely because the median operation obeys the maximum–minimum principle.
In practice, amoeba filters are computed in a spacediscrete setting such that the effective range of spatial frequencies is limited by the sampling theorem. For fixed amoeba radius as in Figure 14, the relevant range of frequencies is determined by the mesh size of the pixel grid. If this mesh size is not below approx. , the higher lobes of the amplification function that make up the difference to selfsnakes with presmoothing do not take effect at all. Translating this to a grid with mesh size , as common in image processing, this means that for amoeba radius up to approx. the frequency response of amoeba median filtering does almost not differ from that of selfsnakes with presmoothing.
Test Case 2: LevelLineAligned Oscillation
To complement the perturbation analysis of gradientaligned oscillations, a second test case is considered in which the perturbation frequency is aligned with the level line direction, . The resulting input signal, compare the schematic representation in Figure 12(c), reads
(39) 
This test case was not presented in [54, 57]. Given that selfsnakes act smoothing along level line direction, it can be expected that this kind of perturbation is dampened by their evolution. This will be confirmed by the analysis, and the corresponding behaviour of the amoeba median filter will be stated.
SelfSnakes Analysis Unlike for the first test case, gradient directions of now vary across the image range, combining constant with . Accordingly, the edgestopping function takes the values
(40) 
and thereby , .
This leads further to
(41)  
(42)  
(43) 
thus after inserting into (35)
(44) 
which confirms by the negative sign of the frequency response factor that the perturbation is smoothed out by the selfsnakes process.
Presmoothing here leads to
(45) 
which in the further course of the calculation only influences higherorder terms, such that (44) is replicated.
Remark on explicit time discretisations. A difference to the first test case to be noted here is that the negative amplification factor does not depend on . This implies a time step size limit for explicit time discretisations of presmoothed selfsnakes: With denoting the highest perturbation frequency that can occur in the discretised image, given by the Nyquist frequency of the grid ( for spatial mesh size ), the amplification factor within a single time step of size must not become or lower, thus must be observed.
Amoeba Filter Analysis.
To determine the response of an amoeba median filter step to the perturbation (39), we consider again Euclidean amoeba metric and . The image graph of (39), compare Section 2.4, is a developable surface. The amoeba structuring element around then is the projection of a bent Euclidean disk affixed to to the image plane, compare Figure 12(c).
The orthogonal projection of the same bent disk not to the image plane but to the unperturbed image graph , compare Figure 12(d), is symmetric w.r.t. the line ; note that the point projects to