Amoeba Techniques for Shape and Texture Analysis

Amoeba Techniques for Shape and Texture Analysis

Abstract

Morphological amoebas are image-adaptive 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 image-dependent 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 Curvature-based PDE Self-snakes 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 structure-preserving denoising or shape simplification [23]. The fundamental building blocks of classical mathematical morphology are non-linear 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 space-variant 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 amoeba-based active contour method was designed [53, 54, 56]. Recently, a combination of edge-weighted 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 multi-channel 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 image-adaptive structuring elements in the space-discrete as well as the space-continuous 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 space-continuous 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 edge-stopping 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 edge-stopping 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 max-plus/min-plus convolution and conjugacy of structure elements.

Section 4 considers the application of amoeba techniques to devise basic algorithms for unsupervised segmentation of grey-value 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 pre-smoothed versions of these, in order to reduce noise sensitivity of filters and to improve numerical stability. This is also the case with self-snakes and active contour PDEs; note that the self-snakes PDE is even ill-posed without such pre-smoothing. Section 5 investigates the effect of pre-smoothing in the self-snakes 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 multi-channel GAC segmentation.

2 Morphological Amoebas

Well-known local image filters such as the mean filter, median filter, morphological dilation or erosion consist of two steps: a sliding-window 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 grey-value 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 grey-value contrast into a modified metric on the image.

2.1 Edge-weighted Neighbourhood Graph

To define morphological amoebas on discrete images, we start by considering edge-weighted graphs based on the image grid.

Definition 1.

Let be a discrete image. Construct an edge-weighted 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 edge-weighted graph is called neighbourhood graph of .

In this definition, neighbourhood can be understood as a 4-neighbourhood, as done in [32], or as an 8-neighbourhood 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 edge-weighted 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 space-continuous versions of amoeba filters; the results proven there allow to interpret amoeba filters as time steps of explicit discretisations for suitable PDEs.

To devise space-continuous versions of amoeba filters, one has to translate first the notion of amoeba metric to the space-continuous setting. Once this is done, the definition of an amoeba as a -ball around a reference point is straightforward.

The amoeba metric for a space-continuous greyvalue image – a real-valued 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.

unit disk on amoeba inimage planeimage graph image plane
Figure 1: Amoeba as projection of the unit disk on the image graph to the image plane.

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.

(a)(b)(c)(d)(e)(f)
Figure 2: Typical shapes of amoebas in the continuous domain for different amoeba metrics. The left group of three, (a–c), shows amoebas on an image with equidistant straight level lines, the right group of three, (d–f), shows amoebas on curved level lines (schematic). Frames (a) and (d) show amoeba metric, (b) and (e) show Euclidean amoeba metric, (c) and (f) show the maximum amoeba metric. Each amoeba is shown with its reference point (bold) and level line through the reference point (dashed).

Figure 2 shows typical amoeba shapes in smooth image regions for the three exemplary amoeba metrics exposed in Section 2.2.

3 Amoeba-Based 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 non-adaptive, sliding-window 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 non-trivial steady states, so-called 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 non-adaptive 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.

abc
Figure 3: Non-adaptive and amoeba median filtering. (a) Original image. – (b) Filtered by 5 iterations of standard median filtering with a discrete disk of radius as structuring element. – (c) Filtered by 5 iterations of amoeba median filtering with Euclidean amoeba metric, , .

PDE Approximation

As noticed already in 1997 by Guichard and Morel [20], the overall robust denoising effect and the characteristic corner-rounding behaviour of standard median filtering resemble the properties of the well-known (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 so-called self-snakes filter [44] allows curvature-based 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 self-snakes 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 self-snakes PDE [44]

(11)

where is a decreasing edge-stopping 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 area-bisecting 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.

(a)(b)
Figure 4: (a) Amoeba with reference point , level line through (dot-dashed) and bisecting level line (dashed), schematic. – (b) Amoeba with curvilinear coordinate system formed by level lines (dashed) and gradient flow lines (solid).

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 sub-intervals 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).

(a)straightlevellineasymmetricamoeba(b)curvedlevellinesymmetricamoeba(c)
Figure 5: (a) Amoeba with straight level line (dot-dashed) through its reference point and further radial lines (dashed) of a polar coordinate system centred at . (b) Area difference in an asymmetric amoeba with straight level lines. The hashed region is enclosed between the right arc of the amoeba contour and the point-mirrored copy of its left arc. – (c) Area difference in a symmetric amoeba with curved level lines. – (b), (c) from [54].

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. Cross-effects 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 right-hand side of the PDE approximated by the amoeba filter, can be derived.

Amoeba Metrics and Edge-Stopping Functions

It remains to specify the relation between amoeba metric and edge-stopping 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 Perona-Malik diffusivity [39] that is also one of the common choices for in the self-snakes equation.

When using the amoeba metric, , the integral in (14) can be numerically evaluated, and one obtains an edge-stopping 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 edge-stopping functions are depicted in Figure 6.

Figure 6: Edge-stopping functions , and associated to , Euclidean and amoeba metrics, respectively. Throughout these metrics, the contrast scale has been set to .

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 amoeba-based 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 max-plus 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 max-plus 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 min-plus 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 max-plus convolution, the right-hand side (20) is the max-plus analogon of a (discretised) integral operator. Herein, acts as the max-plus counterpart of just the same type of integral kernel that appears as point-spread function in space-variant image deconvolution models.

Similarly, amoeba erosion becomes a min-plus integral operator with a min-plus kernel . Generally, conjugate structuring elements in the space-variant 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 self-conjugacy of amoebas.

Figure 8 shows the results of non-adaptive and amoeba dilation and erosion of an example image depicted in Figure 7. Non-adaptive as well as amoeba-based 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.

Figure 7: Grey-scale image ( pixels) used to demonstrate non-adaptive and amoeba-based morphological filters.
abcd
Figure 8: Morphological dilation and erosion, non-adaptive and amoeba-based, of the test image from Figure 7. (a) Non-adaptive morphological dilation with disk of radius as structuring element. – (b) Amoeba dilation with Euclidean amoeba metric, , . – (c) Non-adaptive morphological erosion with structuring element as in (a). – (d) Amoeba erosion with amoeba parameters as in (b).

PDE Approximation

It is a well-known fact that Hamilton-Jacobi PDEs

(24)

describe dilation (“” case) and erosion (“”) of continuous-scale images or level-set functions in the sense that evolution of an initial image by (24) up to time yields the dilation or erosion of with a Euclidean ball-shaped 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 Hamilton-Jacobi-type PDE

(25)

where the “” sign applies for dilation, and “” for erosion.

The proof of this result can be found in [57]; it is based on Proof Strategy I from Section 3.1.2.

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 Hamilton-Jacobi 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 point-symmetric 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 self-conjugacy (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.

abcd
Figure 9: Non-adaptive and amoeba-based morphological closing and opening applied to the test image from Figure 8. (a) Non-adaptive closing with disk-shaped structuring element of radius . – (b) Amoeba closing with Euclidean amoeba metric, , . – (c) Non-adaptive opening with structuring element as in (a). – (d) Amoeba opening with amoeba parameters as in (b).

In Figure 9 exemplary results of non-adaptive and amoeba-based closing and opening of the test image from Figure 7 are shown. Like its non-adaptive counterparts, amoeba-based closing and opening remove small-scale 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 non-adaptive 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 semi-group 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 (non-adaptive) structuring element radius is equivalent to dilating once with radius . Such an additive semi-group 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 semi-group 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 PDE-based filters may be considered in future research.

4 Grey-Scale 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. Intensity-based 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, contour-based 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 state-of-the-art 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 level-set 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 high-contrast locations. The PDE for GAC in level-set representation reads

(31)

Herein, is the evolving level-set function in the plane that represents the actual evolving contour as one of its level sets (by default, the zero-level set), and is the invariable image being segmented. The similarity of (31) to self-snakes (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:

  1. Compute amoeba structuring elements based on the input image .

  2. Initialise the evolving level-set function to represent the initial contour.

  3. 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 space-continuous 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 level-set 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 near-convergence 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 small-scaled 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.

abcd
Figure 10: Amoeba and geodesic active contour segmentation. (a) Detail ( pixels) from an MR slice of a human brain with initial contour enclosing the cerebellum. – (b) Amoeba active contours with Euclidean amoeba metric, , , 10 iterations. – (c) Amoeba active contours with amoeba metric, , , 60 iterations. – (d) Geodesic active contours with Perona-Malik edge-stopping function, , 960 iterations of explicit scheme with time step size . – From [53, 56].

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 short-cut 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 right-hand side of a Hamilton-Jacobi 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 contrast-dependent, was mentioned. This has been done in [9, 31, 35] by modulating the force term in a geodesic active contour model with the same edge-stopping function , such that the entire force term reads as with constant .

The relation between amoeba quantile filters and Hamilton-Jacobi 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].

Quantile bias. Choosing the element with index from the rank order within each amoeba approximates a force term with . In the rotationally symmetric case this term lies between the constant weight of [10] and the -weight from [31].

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 finite-difference scheme.

abcd
Figure 11: Segmentation with initialisation inside the sought object by amoeba and geodesic active contours with dilation force. (a) Detail ( pixels) from an MR slice of human brain with initial contour placed inside the corpus callosum. – (b) Amoeba active contour evolution with Euclidean amoeba metric, , , fixed offset bias , and 20 iterations. – (c) Same as in (b) but 35 iterations. (d) Geodesic active contours with Perona-Malik edge-stopping function, , dilation force (multiplied with the edge-stopping function) and erosion force (independent of the edge-stopping function), explicit scheme with time step size , 18,960,000 iterations. – From [56].

5 Pre-Smoothing in Self-Snakes and Amoeba Filters

The approximation result of Theorem 1 associates iterated amoeba median filtering with the self-snakes equation (11). Unlike (mean) curvature motion (10), self-snakes possess edge-enhancing properties. Rewriting (11) by the product rule, one can state the self-snakes 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 edge-enhancing behaviour. Unfortunately, this term has a shock-filter 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 ill-posedness 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 self-snakes evolution differs significantly if the underlying grid resolution is changed.

A common remedy to this ill-posed behaviour is to use pre-smoothing in the argument of the edge-stopping function, i.e. to replace in (11) or (35) by where is the result of convolving with a Gaussian of standard deviation . Thereby, the ill-posedness of self-snakes is removed, and a stable filtering achieved, at the cost of the additional smoothing-scale 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 pre-smoothing 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 Pre-Smoothing in Amoeba Median Filtering, and Amoeba Radius

First of all, notice that a straightforward translation of the pre-smoothing 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 self-snakes with pre-smoothing, but closely related to it.

At second glance, however, it can be questioned whether the introduction of the smoothing-scale parameter into the amoeba median filter is necessary. Unlike finite-difference schemes for self-snakes, amoeba filtering by construction already involves a very similar smoothing-scale parameter, namely, the amoeba radius . One can conjecture that the positive necessarily used in any amoeba computation could already provide a pre-smoothing 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 self-snakes and amoeba median filter evolutions, see Figure 12(a). From this slope, described by the function , , test cases are derived by adding small single-frequency 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 first-order effects of perturbations still gives a useful intuition about the behaviour of the filters.

Test Case 1: Gradient-Aligned 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)

Self-Snakes Analysis. To determine the response of the self-snakes 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 stair-casing behaviour of self-snakes without pre-smoothing.

Using pre-smoothed in the edge-stopping function argument, one has instead , and therefore

(38)

with the frequency response factor that is globally bounded with its maximum at . Therefore, pre-smoothing 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 space-continuous 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.

(a)(b)(c) (d)
Figure 12: Schematic representation of example functions used in the perturbation analysis, Section 5.2. (a) Graph of unperturbed function , with a Euclidean -disk whose projection to the - plane yields an amoeba. – (b) Graph of a function of type (36) including a gradient-aligned perturbation. – (c) Graph of a function of type (39) including a level-line-aligned perturbation. – (d) Cut in direction through the graph from (c) and the unperturbed graph from (a). The sketch includes further the amoeba around , the corresponding Euclidean disk on and the projection of to which is centred at .
(amoeba filtered), amplification factor (amoeba filtered), amplification factor
Figure 13: Numerical computation results for the amplification of a gradient-aligned perturbation of a linear slope function by one amoeba median filtering step. Top row shows , bottom row . Graphs in left column show unperturbed function , perturbed input function , and filter result ; graphs in right column show perturbations and . Horizontal axes represent , vertical axes represent function values. Computations were carried out on a grid with mesh size .

Figure 14 shows the amplification function for together with its counterpart for one time step of self-snakes with pre-smoothing, 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 self-snakes 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 pre-smoothed self-snakes equation and amoeba median filtering. However, for higher frequencies the amplification factor of pre-smoothed self-snakes 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 pre-smoothed self-snakes 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 space-continuous 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 space-discrete 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 self-snakes with pre-smoothing 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 self-snakes with pre-smoothing.

(amoebas, )(self-snakes, )computed (amoebas, )
Figure 14: Amplification of a gradient-aligned perturbation of a linear slope function by one amoeba median filtering step (theoretical and numerical values) and a corresponding time step of an explicit scheme for self-snakes with pre-smoothing. Adapted and extended from [54].

Test Case 2: Level-Line-Aligned Oscillation

To complement the perturbation analysis of gradient-aligned 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 self-snakes 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.

Self-Snakes Analysis Unlike for the first test case, gradient directions of now vary across the image range, combining constant with . Accordingly, the edge-stopping 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 self-snakes process.

Pre-smoothing here leads to

(45)

which in the further course of the calculation only influences higher-order 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 pre-smoothed self-snakes: 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.
(amoebas, )(self-snakes)computed (amoebas, )
Figure 15: Amplification of a level-line-aligned perturbation of a linear slope function by one amoeba median filtering step (theoretical and numerical values) and a corresponding time step of an explicit scheme for self-snakes (with or without pre-smoothing).

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 . Moreover, the projection from to changes areas only by a factor