2

Literature Survey


 

The problem of measuring the size of a leg ulcer given an image may be cast as the equivalent problem of segmenting the wound object from the image background. The immediate background in this case will be the surrounding skin. It is also common for the edge of the limb or foot to be visible so that the limb itself is bordered on either one or both sides by the background to the entire scene. This chapter begins by assessing the image properties required by various segmentation algorithms for good performance and contrasts these with the properties that are common to wound images. The second part of the chapter considers active contour models in detail, considering the problems associated with the original model due to Kass et al. (1987) and the modifications and alternatives reported by later works. Lastly, the chapter discusses the problem of initialisation of active contour models, which if automated could generalise the algorithms proposed in this thesis to be automatic rather than semi-automatic.

 

2.1 Image Segmentation Methods

The aim of image segmentation is to partition an image into a number of non-overlapping regions that form a complete tessellation of the image plane. A wide range of work has been undertaken to achieve this aim and segmentation has found diverse applications ranging from medical to military uses. It is still a subject of on-going investigation and it cannot be conclusively stated that the segmentation problem has been solved. For the goal of measuring the area of a leg ulcer it is not necessary to obtain a complete segmentation of the image, but to obtain a partial segmentation that discriminates between wounded and non-wounded regions. Surveys of segmentation methods concerning the most common algorithms have been undertaken by Fu and Mui (1981) and Haralick and Shapiro (1985) and most good image processing texts generally provide a broad over-view of the more common algorithms, e.g. Gonzalez and Woods (1992) and Ballard and Brown (1982). The following discussion appraises the suitability of published algorithms under several sub-headings, although there is inevitably some overlap between them.

 

2.1.1 Histogram Analysis

The most basic segmentation procedure that may be carried out is thresholding of a gray level image – see Sahoo et al. 1988. The success of this approach hinges upon whether suitable thresholds exist and whether they may be inferred from the image histogram. Thresholds may be applied globally across the image (static threshold) or may be applied locally so that the threshold varies dynamically across the image. Under controlled conditions, if the surface reflectance of the objects or regions to be segmented are uniform and distinct from the background and the scene is evenly illuminated then the resulting image will contain homogenous regions with high-definition boundaries that generally lead to a bimodal or multi-modal histogram. Various methods exist for determining an appropriate threshold, e.g. Weszka (1978), Kohler (1981) and Kittler and Illingworth (1985). For the more general case, where surface reflectances vary and the illumination is not uniform, the threshold needs to vary across the image. An example of a spatially varying threshold based on local histogram analysis is Chow and Kaneko (1972). This algorithm computes a variable threshold function over the whole image by fitting minimum-error thresholds to gray level histograms of local candidate windows.

A novel two-stage threshold process is used by Wu et al. (1995) for segmenting cell images. The authors state that the images used suffer from low contrast and variable intensity and are unevenly illuminated. The region of interest containing the cell is first roughly located by applying a minimum error threshold (Kittler and Illingworth, 1986) to an image of the local gray level variance. A second threshold due to Otsu (1979) is then applied to the local region only, which separates the cell from its immediate background. Morphological closing and hole-filling operations are subsequently required to remove segmentation artefacts present within the cell region and at its boundary. This appears to be a robust approach to image segmentation that shares some problems common to segmenting leg ulcers. However, the local gray level variance of a leg ulcer image is dependent upon the local curvature of the limb, specular reflections and slough. Figure 2.1 shows example gray level histograms produced from two leg ulcer images.

 

(a)

(b)

Figure 2.1 Example gray level histograms for (a) image 1.1 in Plate 1. The surrounding skin tissue in this image appears relatively unaffected and is thus quite homogenous. An area of interest was selected that eliminated the background. (b) An example of a wound image (image 2.1 in Plate 2) which lacks homogeneity of gray scale.

 

Figure 2.1(a) shows that a shallow valley point exists in the histogram of the image gray level, implying that a threshold could be quite successful in discriminating the wound from the surrounding skin. The images shown in Plates 1 and 2 are twenty examples extracted from a wound library that has been compiled from hospital visits, a proportion of which have been captured using the MAVIS instrument (Plassmann, 1992). Many wound images show substantial changes in gray level as a result of practical difficulties that arise in obtaining well-illuminated images. Difficulty in evenly illuminating highly curved parts of the leg is a common cause of this problem. Clearly, this problem is exacerbated as the physical size of the wound increases. Further difficulties exist in accurately locating the wound boundary when epithelialisation tissue is present at the wound edge, examples of which are images 1.3, 1.4 and 1.5 in Plate 1 and image 2.7 in Plate 2.

 

2.1.2 Region Growing and Split-and-Merge Algorithms

Region growing algorithms start with a number of seed pixels or seed regions and grow these regions by adding to a region previously unassigned neighbouring pixels that correspond to some similarity criterion for that region (Gonzalez and Woods, 1992). The rule employed for accepting/rejecting a candidate pixel for merging with the neighbouring region may be defined as either on the basis of a pixel and its immediate neighbours or on its relationship to the mean of a region (using either Euclidean or Mahalanobis distance as the discriminant). Growing algorithms that link neighbouring pixels by analysing contrast use edge-strength to measure the compatibility between neighbouring pixels (edge strength) and thus require a completely closed boundary of high-edge strength pixels, otherwise neighbouring regions will be merged. Growing algorithms that measure similarity by comparing a pixel with some measure of the whole region (e.g. mean value) are more robust but impose a stronger homogeneity criterion upon the regions being grown. Recently, Adams and Bischof (1994) presented a region-growing algorithm that is robust and free of tuning parameters. However, they state that the algorithm requires suitable pre-processing to remove background variation due to uneven illumination. Another recent development of the region-growing paradigm is the ‘statistical snake’ due to Ivins and Porril (1994). This algorithm imposes a smoothness constraint on the boundary of a region. Snakes (active contour models) are discussed in §2.2.

Split-and-merge algorithms have been reported, among others, by Klinger (1973) and Horowitz and Pavlidis (1976). This method proceeds to successively divide an image into smaller non-overlapping regions if some similarity criterion is not met, otherwise no split of that region is made. The end result of the splitting is an over-segmented image. A merging procedure is then applied to merge neighbouring regions under the same homogeneity predicate that was used for splitting. Homogeneity tests may be based on simple image statistics, e.g. the range of intensity values in a region, or upon more robust statistical similarity measures such as the Kolmogorov-Smirnov test (Muerle and Allen, 1968). Chou et al. (1992) use an adaptive split-and-merge method to segment magnetic resonance (MR) and x-ray computer-tomography (CT) images for calculation of volume and reconstruction of internal organs. They use a likelihood ratio test based on a Gaussian distribution assumption to check for region homogeneity. The lack of homogeneity across the wound in many images of leg ulcers and the variable contrast make split-and-merge schemes unsuitable for measuring the area of many wounds (see Plates 1 and 2 for examples of leg ulcer images).

 

2.1.3 Colour-Space Analysis and Segmentation

Recursive region-splitting algorithms using colour-based discriminants have been described by Ohlander et al. (1978) and Ohta et al. (1980). Starting with the whole image as one region, these algorithms apply a set of thresholds derived from analysing histograms of individual colour features. Once a region has been segmented, the algorithm is recursively applied to each new region. Ohta et al. examined a range of images and determined a set of three colour features which had near-optimal discriminant power. The primary discriminant in the cases they analysed was found to be equal to the gray level or intensity component. Arnqvist et al. (1988) perform colour analysis of wounds using Red/Yellow/Black classification, corresponding to granulation tissue, slough and necrotic material respectively. Colour images of leg ulcers are automatically analysed by Beriss and Sangwine (1997) using a technique based on 3-D histogram clustering of the colour values. Clustering of the histogram is achieved by erosion and dilation operations that break the connections between the clusters corresponding to different tissue classes. The result of segmenting one wound image into regions of granulation tissue, slough and surrounding skin is presented.

Umbaugh et al. (1993) present two colour segmentation algorithms, which they use to segment images of skin tumours. The first algorithm represents the chromatic component of an image by two angles representing the orientation and elevation of a colour vector in the RGB colour space. This is equivalent to normalising the colour vector using the L2-norm (Healey, 1992). The second algorithm recursively applies a principal component transform (PCT, also known as an eigenvector transform or the discrete Karhunen-Loeve transform (KLT)) to the colour data followed by a split of the colour space at the median point along the principal axis. This process is repeated on each subspace so formed until the colour space is divided into the desired number of colours. A segmentation based on global colour properties as produced by an RGB colour histogram assumed that the colour properties of the wound are homogenous under the combination of illuminants used and level shifts and gamma-transformations applied by the camera and image-capture equipment.

Schmid and Fischer (1997) segment pigmented skin lesions using a principal component transform of RGB colour images of digitised epiluminescense microscopy slides. The 2-D histogram of the first and second principal axes of the data is smoothed by a morphological process before being clustered into three classes of tissue using the fuzzy k-means algorithm. One example of the smoothed 2-D histogram is presented showing good cluster separation. The segmentation result is affected by hair and the inherent texture of the lesion and so a median filter is applied before image segmentation and the subsequent boundaries detected are smoothed and small regions removed by morphological opening and closing.

 

2.1.4 Pixel Classification for Segmentation

Image segmentation can be achieved by grouping together neighbouring pixels that are similarly classified. The classification may be based on trained parameters (supervised classification) or based on clustering techniques, e.g. the k-means and ISODATA algorithms (Duda and Hart, 1973). The automated colour-image analysis mentioned above is an example of unsupervised classification or clustering. Classification schemes can be potentially improved by using contextual information, basic examples of which are relaxation labelling (James, 1987) and Markov random field methods (Therrien, 1989). Däschlein et al. (1994) used an artificial neural network (ANN) to perform the classification task for segmenting CT images. The features assigned to each pixel were based on five basic measures of a local co-occurrence matrix. A contextual classification of medical images is reported by Ossen et al. (1994). They use an ANN to estimate the a-posteriori Bayesian probabilities for a minimum-error classification to aid in the diagnosis of Graves’ ophthalmopathy. Jackson et al. (1993) compared the performance of a relaxation-labelling algorithm with a minimum Mahalanobis distance classifier to segment MR and CT images. They report that the relaxation scheme is more consistent in labelling regions of gray and white matter and cerebro-spinal fluid (CSF). Whilst the visual appearance of the results is impressive it is necessary to consider that the regions of the ‘gray level’ image corresponding to particle densities are somewhat distinct and homogenous. This is in contrast to the conditions found in many images of leg-ulcers.

 

2.1.5 Gradient Operators, Edge and Zero-Crossing Detectors

Figure 2.2 Example of edges detected by applying a threshold to a non-maximally suppressed gradient map produced using a derivative of Gaussian operator with filter parameter s =8 pixels.

Canny (1986) proposed the derivative-of-Gaussian filter as a near-optimal filter with respect to three edge-finding criteria: (a) good localisation of the edge, (b) one response to one edge and (c) high probability of detecting true edge points and low probability of falsely detecting non-edge points. Deriche (1987) implements a filter with an impulse response similar to that of the derivative of Gaussian, but which lends itself to direct implementation as a recursive filter. This requires fewer coefficients than a transversal or non-recursive filter that has the same impulse response, yielding faster convolution, especially when the width of the filter becomes large. The benefit of these filters is that they impose a smoothing upon the image and as such are band-pass filters rather than simple differentiators and thus they are less sensitive to noise. Smoothing is controlled by a scale parameter which may be tuned for each application. Edges are extracted from the gradient image by choosing pixels with a gradient magnitude that is a maximum in the local direction of the gradient.

Marr and Hildreth (1980) have proposed finding zero crossings produced by the Laplacian-of-Gaussian (LoG) operator for the detection of edges. This operator is orientation independent and accurately locates zero crossings, provided that the smoothed gray level variation along, but not necessarily near to, a line of zero-crossings is linear. Haralick (1982) proposed fitting bi-quadratic facet models to an image and applied a zero crossing detector to the second derivative of the function taken in the direction of maximum gradient. This is a more computationally expensive method than the LoG operator proposed by Marr and Hildreth.

In relation to finding the edges of wounds, the application of such derivative operators generally fails to produce closed boundaries – this applies both to maximum gradient detectors and zero crossing detectors. Although it is possible to close small gaps by extending the ends of detected edge chains (Niblack, 1987) it is quite common for edge segments to curve into or away from the wound boundary, even at high levels of smoothing, so that one end of the detected edge segment is positioned some distance away from the wound boundary. In Figure 2.2 the detected edge contour deviates outside of the wound at the lower left side of the wound. This image is shown in colour in Plate 2 (image 2.7). However, for a broad range of leg ulcers, the information provided by a gradient operator provides a strong cue for detecting the edge of a wound. Clearly, such operators will detect all edges and not just those belonging to a wound.

Edge focusing algorithms (Bergholm, 1987) use smoothed gradient operators and detect and track edges throughout a scale hierarchy. There is a much higher likelihood of locating the boundary once the gradient image has been smoothed to a relatively high degree. Typically, the smoothing parameter is reduced by a small amount so that the boundary moves by only a small amount (e.g. 1 pixel), allowing the algorithm to accurately track the movement of the contour detected at the previous scale. Focusing (de-blurring) allows the contour to converge towards its true location in the original unblurred image, since the smoothing displaces and smoothes apparent contours in the image. Denton et al. (1995) use an edge focusing algorithm to detect skin lesion boundaries. Accurate location of the boundary is important so that automated diagnosis can subsequently be made using colour, texture and shape measures to classify the lesion, which may be cancerous. The approximate size of the lesion is assessed by an automated bi-level thresholding operation (Kittler et al., 1984) after pre-processing by a median filter. The estimated lesion size is used to tune the initial scale-constant of the LoG filter used for detecting lesion edges. The success of edge focusing relies upon detecting most of the required boundary at high scales. If any gaps in the boundary were to be left unlinked then neighbouring ends of the edge chains could migrate as the image was de-blurred. Such edge chains occur infrequently in wound images thus limiting their applicability within the scope of leg ulcer measurement.

Perona and Malik (1990) propose a multi-scale or scale space representation of an image based on anisotropic diffusion. This models the image smoothing as the diffusion of heat energy in a physical system. The heat conduction coefficient is allowed to vary across the image plane and is made dependent upon the image gradient. This effectively leads to a spatially adaptive smoothing which tends to preserve the location of edges throughout the scale hierarchy. The idea behind this method is that it avoids unconstrained smoothing of the image across edges, which causes the edges to move and thus prevents them from being accurately located at high scales. The problem formulation leads to an iterative minimisation algorithm and thus the algorithm imposes an increased computational cost.

 

2.1.6 Watershed Segmentation

Figure 2.3 Over-segmentation of an image using watershed algorithm (produced with Visilog image processing package). The gradient was computed using a derivative of Gaussian operator with filter parameter s =8 pixels.

The watershed algorithm may be applied to a gradient magnitude image for the purposes of image segmentation (Haris et al., 1998 and Gauch 1999). Considering the gradient magnitude at each pixel as the height of a surface in 3-D, regions are formed by simulating a flooding of the image that begins at local minima of the gradient image function. The boundaries of each region stop advancing when neighbouring flooding regions meet and the boundaries thus correspond to ridge-paths of the gradient magnitude. The principal advantage of the watershed segmentation scheme over other edge-finding schemes is that it generates closed boundaries. The regions defined by the closed boundaries represent an over-segmentation of the image, since there will usually be many such regions corresponding to each object. If the gradients are computed at successively higher scales, the number of local minima (flood basins) in the gradient magnitude image will decrease. Watershed segmentation can be applied to wound images and will form a sensible over-segmentation of many of the moderate and well-defined images. Figure 2.3 shows an example of applying the watershed segmentation algorithm to a derivative of Gaussian gradient magnitude image. Although the watershed algorithm can produce boundary segments that correspond to each part of the true wound boundary (approximately), the problem of merging wound regions into one region is confounded by the fact that very often two regions, one on either side of the boundary, are similar in their gray level properties and are very likely to be merged. A possible improvement in comparing regions for merging is to make use of second order image statistics in the form of a co-occurrence matrix, as proposed by Chen and Pavlidis (1979).

 

2.1.7 Radial Search Algorithms

Radial search algorithms simplify the task of edge detection by constraining the search for an edge point to lie on a line emanating radially from a central point within the object. The method also simplifies the detection of erroneously marked points by allowing a straightforward comparison between edge points marked on neighbouring radial lines. Golston et al. (1990) used a radial search algorithm to detect skin tumour boundaries having considered that standard methods were only suitable for a small sub-set of tumour images. A heuristic measure for edge-worthiness along each radial search line requires the luminance to increase by a specified minimum amount for a specified minimum sustained number of pixels. The authors reported good results for 17 out of 20 tumour images tested. Radial search algorithms are only applicable to delineating objects that are mainly convex. Ruiz and Fairhurst (1995) detected cardiac boundaries in 2-D echocardiographic images by applying a LoG filter and zero-crossing detector to the image prior to radial sampling. O’Gorman et al. (1983) use a radial search to segment cell nuclei. On each radial line a sigmoidal-shaped monotonically increasing edge model is used to select an appropriate gray level threshold. Although such algorithms have been applied to difficult segmentation tasks with low and variable boundary contrast, leg ulcers give rise to many potential edge-point candidates within the wound itself.

 

2.1.8 Boundary Following Algorithms

Boundary following algorithms have been used to find the edges of cells in microscopy images. An algorithm due to Jarkrans et al. (1980) segments cell nuclei and cytoplasm by searching for gradient maxima along 8 radial lines emanating from a manually selected seed point. The points detected by the radial search are joined by a tracking procedure that attempts to follow a path of maximal gradient given by a chain-code of the gradient direction. Not all attempts to track cell-boundaries were successful in producing fully closed boundaries. Graph-searching techniques for boundary detection based on a heuristic graph-searching algorithm by Nilsson (1971) have been reported by Martelli (1972) and Lester et al. (1978).

The principal disadvantage of boundary-following algorithms is that enough boundary information needs to be available in the image to form an acceptable path between start and end search nodes. Expanding all possible search nodes between two points does not guarantee that a path will be found linking the two. As already stated, the gradients determined for leg ulcer boundaries form incomplete paths, are ambiguous and can trace their way into the wound rather than following the actual boundary.

 

2.1.9 Summary

This section has discussed various methods of image segmentation and some of their applications. Many of the methods could be used to segment some of the images that are more evenly illuminated, free from specular reflection and more homogenous with regard to their gray level properties. However, such wounds are fairly uncommon among the wound images observed and in general any segmentation algorithm will produce progressively less satisfactory results as the ideal conditions for an image degrade. The vast variability in the quality of wound images with respect to these idealised characteristics generally prevents any one algorithm from being universally applied. If an algorithm is to be successfully applied to measuring the area of leg ulcers then it has to be capable of producing sensible results on the majority of wound images. Clearly, the accuracy of such an algorithm will be affected by the qualities of the wound image. However, if the algorithm can still produce a reasonable measurement result when the image quality is poor then this represents a best solution under the circumstances. It is this variability which precludes the general use of any algorithm which relies upon some quasi-homogeneous property to distinguish the wound from the surrounding tissue. Most of the images of leg ulcers in the library show variation in boundary contrast.

 

2.2 Active Contour Models

Active contour models (also known as ‘Snakes’) were introduced by Kass et al. (1987) as a novel solution to the low-level imaging task of finding salient contours (e.g. edges and lines) in digitised images. An interesting and powerful property of an active contour model is its ability to find subjective contours and interpolate across gaps in edge chains. An active contour model represents an object boundary or some other salient one dimensional image feature as a parametric curve that is allowed to deform from some arbitrary initial position towards the desired final contour. The problem of finding this final contour is cast as an energy minimisation problem with the intention that it yields a local minimum of the associated energy functional. The energy functional is thus based upon the salient features of the image under inspection so that when the contour is aligned with such a feature, the energy of that contour corresponds to a local energy minimum. The minimisation of the energy functional gives rise to spatial forces that act upon the contour to deform it by translating, stretching and bending it as necessary. It is necessary to include in the energy functional terms that impose smoothing constraints upon the contour since the image itself is corrupted by noise from various sources. The original model incorporates two internal energy terms acting to constrain the solution by imposing conditions related to contour smoothness and regularity. These properties are discussed in §2.2.3. Implementation of an active contour model involves several basic considerations that are explored. All aspects of the original model have been subject to modifications and extensions by later researchers that have served to address perceived weaknesses of the complete implementation.

 

2.2.1 Modelling of the Contour

An active contour model algorithm comprises the fitting of a deformable model of either a closed-boundary or open-boundary-segment to available image data, e.g. the outline of an object and in the case of the present work, the edge of a wound. The contour is defined in the (x,y) plane of an image as a parametric curve, v(s)=(x(s),y(s)), where s is a parameter advancing around/along the contour and is related to arc length. This is a continuous-domain representation of the contour. Subsequently it will be necessary to represent the curve in a discretised form where v(s) is either (a) uniformly point-sampled in the domain of s (finite difference discretisation) or (b) represented by a sequence of piecewise polynomial segments (finite element discretisation). It is desirable in this instance that the curve possesses continuous first or second order derivatives across the boundaries of neighbouring elements – referred to as C1 and C2 continuity respectively. The discretisation step is considered in §2.2.6.

 

2.2.2 Formulation of the Energy Functional

Having specified the contour as v(s), the deformable model is defined as a sum of energy terms in the continuous spatial domain, where the energy terms may be categorised as (a) internal, i.e. functions of the contour v(s) itself, or (b) external, generally functions derived from the image under inspection. The original model defines a third category, constraint energy, consisting of higher-level factors added to provide more control over the contour than the internal energy terms, although such factors have not been used in this work. The continuous-domain energy functional for the Kass snake with one external energy term and without constraint terms is given by:

(2.1)

where a and b are coefficients governing the elasticity and stiffness of the contour v(s) respectively and P(v) is an external energy (potential) function sampled along the contour. The convention is adopted of denoting a derivative or partial derivative of a function by appending the appropriate subscript, e.g. vss=d2v(s)/ds2.

 

2.2.3 Minimisation of the Energy Functional

Generally the potential function P(x,y) will give rise to many local minima of (2.1). Any contour v*(s) that corresponds to a local minimum of (2.1) is a solution to the contour-finding problem. Associated with such a functional is a set of governing differential equations, termed the Euler-Lagrange equations, the solution of which yields the energy-minimising contour. Farlow (1982) gives a good introduction to the task of deriving such equations. Minimisation of the energy model of (2.1) gives rise to the following pair of differential equations, the terms of which equate to the various internal and external forces that are balanced when (2.1) is locally minimised:

(2.2)

where v*(s) is a local minimiser of (2.1).

Substituting an arbitrary contour v(s) for the minimising contour, v*(s) into (2.2) generally yields a non-zero residual term. It is this residual function which is responsible for deforming a contour v(s) towards the minimising contour v*(s).

It will help to elucidate the nature of the component energy terms and respective forces by way of some elementary analyses.

 

2.2.3.1 Elastic Energy and Forces

The elastic energy component of the active contour energy functional (2.1) and its resultant forces taken from the Euler-Lagrange equations (2.2) are given respectively by:

(2.3)

(2.4)

Figure 2.4 Elastic forces arising from a closed contour under tension.

The energy term (2.3) is minimised (to zero) when the first derivative of the contour is identically zero everywhere, signifying that the contour of minimum elastic energy is one that has collapsed to a point. One may wish to specify end conditions for a contour, which can be of open or closed (cyclic) form. Consider the case that arises when the contour is an open curve with fixed end points (hard constraints upon (2.1)): clearly the constraints prevent the contour from collapsing to a point and the action of the elastic forces is to straighten the curve. Such action requires forces that are generally directed into the centre of curvature. When the contour is closed the action of the elastic forces is unchanged, although ultimately the contour is deformed to a completely convex object and the forces become directed towards the inside of the object. Figure 2.4 shows a parametrically-defined closed contour upon which are superimposed the force vectors resulting from the second derivative (elastic) term of (2.3). Since the force vectors are always oriented in the direction of the centre of curvature at each point on the contour, at a concavity the forces are oriented outwards from the inside of the object. This causes the contour at these points to expand until the concavity is eliminated. The force vectors in Figure 2.4 are not normal to the boundary however, indicating that the movement of the boundary into the direction of curvature is only one mode of movement. The secondary action of elastic forces is shrink the curve – this action is tangential to the contour and in the case of the open curve is the force component responsible for collapsing it to a point in the absense of competing or constraining forces. The decomposition of elastic force into its normal and tangential components is considered in (§3.2).

 

2.2.3.2 Bending Energy and Forces

The bending energy term and corresponding force term of a contour v(s) are given respectively by the following:

(2.5)

(2.6)

Figure 2.5 Two minimum energy contours constrained to pass through a set of three arbitrary points (marked by circles).

Sharp corners or points of high curvature and fine contour detail are characterised by high frequencies. In comparison with the elastic energy term, the bending energy term is more sensitive to these higher-frequency components because it is based on a higher-order derivative. Figure 2.5 shows a comparison between two open contours that are constrained to pass through 3 points. Under this constraint, the contours are in equilibrium (energy minimised) for (a) an elastic energy term only and (b) bending energy term only. At the central point in Figure 2.5 the contour’s first derivative is discontinuous and its second derivative is thus infinite. The distribution of the bending energy of such a contour would be infinite and also completely concentrated at this one point. The lack of high curvature and sharp corners implies low bending energy and the curve is smooth and C1 continuous everywhere. However, as with the elastic energy term, (2.5) is scale-dependent and the bending energy of a contour can always be reduced, not just by attempting to distribute the curvature evenly along the contour, but also by a contraction of the whole contour. Without opposing forces, the bending energy term will cause a contour to naturally degenerate to a straight line. This form of (2.4) in fact dictates that the x(s) and y(s) components of v(s) both be linear, rather than simply being a linearly related arbitrary form. When the contour is closed the natural degeneracy of the contour under the minimisation of (2.4) is a point since the contour is constrained from becoming linear. This form of degeneracy is a manifestation of the scale-dependent nature of the energy term.

 

2.2.4 Weaknesses of the Kass Model

The original active contour model suffers from the following weaknesses:

Several authors have proposed various methods of overcoming these difficulties and such modifications and improvements are considered in the following section. However, it must be stated that none of the modifications succeed in overcoming the lack of globality of the solution, requiring good initialisation in order to produce a satisfactory solution.

 

2.2.5 Solution of Euler-Lagrange Equations

When the active contour model is not in a stable equilibrium state (corresponding to an energy minimum) the Euler-Lagrange equations of (2.2) result in a set of unbalanced forces so that the right-hand side (RHS) will not be zero. The solution may be sought by setting the RHS to a negative time-varying residual force term. This is similar in operation to dissipating some of the energy of the system via a frictional or resistive component and is the method that was used by Kass et al. (1987). A physical interpretation of the system has been provided by Terzopoulos (1987) that includes an inertial (mass) term allowing a full physical interpretation of the model including the dynamics of the solution. The inertial component allows the model to overcome the attractive forces acting in the shallow energy wells of undesirable and noisy ‘solutions’. In the continuous time and space domains the full inertial model for system motion may be defined by:

(2.7)

where m is the mass density of the contour and g is a damping constant.

This system is in a stable equilibrium state when the time-derivative dependent terms of the RHS decay to zero so that (2.7) corresponds to (2.2). Equation (2.7) may be suitably discretised in both temporal and spatial domains and subsequently solved by iteration through discrete time. Setting m to zero in (2.7) yields the original solution form that corresponds to a mass-less contour which thus lacks the ability to overcome weak local minima. Now, any physical system which stores energy in more than one medium and which is not at rest has the potential to undergo oscillations, where energy is cyclically transferred from one storage medium to another. A simple example is a swinging pendulum where the pendulum’s energy is cyclic exchange of potential and kinetic energy. Setting g to zero in (2.7) corresponds to a non-dissipative or Hamiltonian system and in this case any energy ‘lost’ by virtue of the action of the forces on the left-hand side (LHS) is transferred to the inertial term. None of this energy is dissipated or irrecoverably lost from the system, hence the total energy of the closed system (no energy transfer either in or out of the system) is a constant. Thus, even if the contour v corresponded exactly to the desired energy minimum – both sides of (2.7) being zero – the system would not be at rest. For the system to stabilise and converge to a unique energy minimum it must lose energy by dissipating it via the frictional term. A physical analogue of the original evolutionary equation may be made by considering the mass-density parameter m to be very small. Note that letting m=0 implies that g¹0, otherwise the physical analogue breaks down. The Rayleigh dissipation constant g plays an important role in the stability of the system. If g>>m the system is clearly over-damped and in theory should converge to an energy minimum with no overshoot or oscillation. However, the complex topology of the external potential function P(x,y) dictates that the equations must be solved numerically by discrete time stepping. Transforming the continuous-time system of (2.7) into a discrete time system introduces the possibility of instability.

 

2.2.6 Spatial Discretisation

The differential equation of (2.7) may be spatially discretised by representing the contour v(s) and its derivatives using either finite differences or finite elements. The matrix coefficients arising from finite differencing of (2.7) are given by Kass et al. (1987). In this application it is proposed that finite element discretisation be used. The application of the finite element (FE) method for solving differential equations to the problem of discretising and solving the continuous-domain energy formulation of active contour models has been reported by Karaolani et al. (1990) and Cohen and Cohen (1993). Karaolani et al. compare the use of parabolic and hermite cubic polynomials for spatially discretising an active contour model. The derivatives of a parabolic curve vanish after the third derivative making them unsuitable for active contour models requiring a bending energy component. It is possible, by use of the Rayleigh-Ritz method or by judicious use of integration by parts, to reduce to two the number of derivatives required by a finite element implementation that contains a bending energy term. However, cubic elements provide a method for applying C1 continuity to a curve, which is a desirable aspect when the curve is required to be smooth. They also allow more flexibility in their application.

The finite element approach has the following advantages over finite difference schemes:

 

2.2.6.1 Finite Element Shape Functions

The most convenient form of interpolating curve for this application is that of piecewise-polynomial curves, provided they possess at least C1-continuity to produce smooth boundaries, are twice differentiable and that changes to any one of the control vertices only affects the piece-wise curve locally. Both cubic Hermite and cubic B-Spline curves possess these qualities (Bartels et al., 1987) and are readily used in the formulation of a finite element active contour model, preferably with some adjustment in the case of Hermite curves. Discretising the parametric curve v(s) into finite elements leads to an expression for each element given by:

(2.8)

where N(s) is the set of shape or basis functions defining the form of the curve and Ve is a two-column matrix of control vertices.

The use of both Hermite cubic and cubic B-Spline elements is considered next. Definitions of the respective shape functions are given in Appendix A.

 

Cubic Hermite Basis

Interpolation with a cubic Hermite basis requires specification of both end-points of the curve and both end-derivatives of the curve. This means that the matrix equation for the curve is not directly expressed in the form required for application with the finite-element method. Such curves have great application in the field of computer graphics and are particularly useful for producing inter-frame interpolation for animated sequences (Bartels et al., 1987). In such applications the need to specify the end-point derivatives of the curve segments is a powerful tool for applying tension to the interpolating curve and also allows a bias to be added so that derivatives may be skewed. Furthermore, one may desire to break the implicit C1-continuity at the knot points between piecewise curve segments, so producing a corner, by allowing the two derivative values at this point to differ. Thus, in general, there is good reason to retain control over the end-point derivatives. Hermite curves have the property of interpolating their control vertices.

 

Cubic B-Spline Basis

The set of cubic B-Splines has the property of producing C2-continuity between curve segments (elements). It should be noted that cubic B-Splines do not interpolate the control vertices. To fit a unit-spaced B-Spline that interpolates a set of element end-points Pe one must choose control vertices Ve such that:

(2.9)

Thus, given an approximate wound boundary suitably represented by a vector of element end-points Pe one must obtain the control vertices Ve that interpolate the end points and perform iterations of the active contour model until convergence is attained. Subsequently, the final element end-points may be calculated from the Ve that represent the converged contour. The matrix of equations defined by (2.9) is tri-diagonal and diagonally dominant, allowing rapid solution.

 

Comparison of Hermites and B-Splines

Figure 2.6 shows a comparative plot of the cubic Hermite and cubic B-Spline shape function sets, produced by plotting the respective shape equations listed in Appendix A. The use of cubic B-Splines gives rise to smoother curves (by virtue of higher-order continuity between elements), but such a curve does not interpolate its control vertices, as stated above. The higher order continuity between the elements of a cubic B-Spline curve tends to produce smoother curves than those defined by cubic Hermite elements.

Figure 2.6 Comparison of cubic B-Spline and default-slope cubic Hermite basis (shape) functions

 

2.2.6.2 Stiffness Matrix Equation

The discretisation of the static problem given equivalently by (2.1) and (2.2) yields a matrix equation. At a unique and stable equilibrium, the complete discretised system of equations may be written in matrix form as:

(2.10)

where V* represents the control vertices defining the curve of specified form which locally minimises (2.1).

The derivation of this form of equation from (2.1) is given in Appendix A for both cubic hermite and cubic B-Spline elements. The size of the K matrix, its sparseness and the precise values of its coefficients are dependent upon the number of points or elements (and their form) used to represent v(s) discretely. The cubic hermite curve gives rise to coefficients of significantly larger magnitude than those of the B-Spline, particularly in regard to the stiffness terms. Leymarie and Levine (1993) show that the conditioning of the matrix equations is partially dependent upon the size of the eigenvectors of K which in turn are affected by the magnitude of the coefficients. However, this is unlikely to be a matter of concern unless high values of b are used and the equations are solved with large time steps.

 

2.2.7 Time Discretisation

A direct single-step analytical solution of (2.10) for V* is not feasible since the external forces which give rise to the F(V) term are non-linear. Therefore the aforementioned discrete time stepping method of solution is necessary. When the system of equations is not in equilibrium, the RHS of (2.10) is non-zero and thus (2.10) may be generalised to a time-dependent version given by:

(2.11)

Rearranging (2.11) for V(n+1)T yields an iterative equation:

(2.12)

Equations (2.11) and (2.12) are essentially the equations used by Kass et al. (1987). The internal forces are linear and specified wholly by the stiffness matrix K allowing implicit or backward Euler steps to be taken with respect to these forces. Explicit or forward Euler steps are taken with respect to the external energy since its form is generally non-linear and cannot therefore be easily expressed by an analytical expression at time t=(n+1)T. The external forces are thus treated as remaining constant during each time-step. This mixture of implicit and explicit stepping is loosely termed ‘semi-implicit’.

The solution of (2.12) at each iteration may be implemented in several ways. Kass et al. (1987) and Cohen (1991) use LU decompositions, while Leymarie and Levine (1993) suggest using either Cholesky decomposition or a penta-diagonal algorithm due to Benson and Evans (1977). Cholesky decomposition is applicable when the matrix to be ‘inverted’ is both symmetric and positive definite (all positive eigenvalues). Press et al. (1992) state that "Cholesky decomposition is about a factor of two faster than alternative methods for solving linear equations." Further considerations for speed improvement focus upon exploiting the sparseness of the stiffness matrix – see Benson and Evans (1977) – and upon how often the matrix changes (e.g. by refining T) and thus how often it requires re-inverting. If the matrix remains constant its inverse may be calculated from any of the above decomposition methods and stored so that subsequent steps require only matrix addition and multiplication. In the present case however, such considerations need not be a matter of great concern since (a) a wound does not require a large number of elements to define its boundary and (b) execution time of a wound-measurement algorithm is not too critical. The MAVIS instrument produces results within five minutes of the wound image being taken (Plassmann, 1992).

 

2.2.8 Stability of the Iterative Solution Method

Numerical solution of equations by implicit stepping methods is only guaranteed to be stable when a system is linear. However, the external force term in (2.12) gives rise to possible numerical instability. If the external forces dominate the internal contour smoothing forces in (2.12) and too large a time-step is taken, the contour could easily be moved out of the domain of influence of a local edge-profile and become influenced to converge to another edge that is not part of the desired solution. Cohen (1991) addresses this situation by normalising the external forces to have unit magnitude and adds a scaling factor, k. The external forces then become ktFt/|Ft| where t=T/g. A second problem Cohen considers occurs when a contour attempts to stabilise in alignment with a path of maximum gradient. Since the external force function is defined on a finite grid, the forces are suddenly reversed as the contour moves under the action of those forces from one side of the edge to the other. Figure 2.7 shows how this can result in a limit cycle oscillation:

(a)

(b)

Figure 2.7 (a) Production of limit cycles from non-interpolated external forces (b) Convergence of steps due to use of linear interpolation.

 

However, this does not complete the analysis since the smoothing applied by (I+tK)-1 can itself cause the contour to move by a large amount when the internal forces dominate the external forces – this is independent of the choice of t. Equation (2.13) analyses the problem by expressing the displacement of the control vertices under the action of (2.12):

 

(2.13)

The term in the square brackets is easily recognisable as the residual force term from setting (2.10) to be non-zero, with the external forces modified as suggested by Cohen. Clearly, both parts of the residual at time t are smoothed by the inverse-matrix term and scaled by the time step to produce the displacement of the control vertices controlling v(s). Clearly, when a and b are large enough the internal forces dominate the equation – but there is no fixed limit over exactly how far the time-step taken will cause vertices and hence the contour to move. The displacement of the control vertices given by (2.13) when the internal forces are dominant is a function of Vt itself, a and b as well as t. Therefore, one cannot ensure that a contour will definitely behave in a steepest descent manner by adjusting the time step alone.

 

2.2.9 Direct Energy Minimisation Approaches

Several authors have proposed iterative methods of directly minimising energy functionals of similar form to (2.1). Amini et al. (1991) employed a Dynamic Programming technique (Ballard and Brown, 1982). This involves evaluating the contour energies produced by local changes to the nodes of the contour and is numerically stable and guaranteed to converge to a unique solution within a local minimum of the energy functional. This approach allows the setting of hard constraints and Amini cites the example of constraining the spacing between contour nodes to be greater than a certain specified distance. However, the dynamic programming method of solution has two main disadvantages: (a) It is relatively slow, executing in O(nm3) time where n is the number of contour sample points (‘snaxels’) and m is the number of trial offset patterns associated with each snaxel; (b) the energy functional must be capable of being decomposed into locally dependent terms so that the energy associated with each snaxel is affected by only a few neighbouring snaxels.

To overcome these disadvantages, Williams and Shah (1992) have proposed a ‘greedy’ algorithm, which considers the energy of a single snaxel at each step and moves it to the neighbouring position that minimises the associated energy term. The greedy algorithm executes in O(nm) time, sacrificing the ability to converge completely to a stable minimum. Lai and Chin (1993) propose a modified greedy algorithm that is guaranteed to converge. This is achieved by minimising of the sum of all the local energy terms that are affected by movement of any single snaxel. The original active contour model uses a semi-implicit scheme to promote stability under all settings of the step size parameter, given some upper limit on the size of the external forces. The trade-off for efficient convergence requires a balance between the size of step a contour may take at any one iteration and the need to control the likelihood of the contour moving out of the neighbourhood of one local minimum of the energy functional and into an adjacent minimum.

 

2.2.10 Element Collapse Problem

One of the problems encountered with active contour models listed in §2.2.5 concerns the collapse of finite elements and the clustering of finite difference snaxels at points of high gradient. The problem has been noted and solutions proposed by both Amini et al. (1991) and Ranganath (1992). Amini’s dynamic programming solution allows one to specify a minimum spacing between contour nodes and Ranganath’s algorithm re-samples the implied contour after every iteration thereby enforcing uniformly spaced snaxels. The problem is caused by the gradient magnitude varying along the contour so that even when the gradient of the external energy term normal to the contour is zero everywhere along the contour, there is a non-zero gradient component tangential to the contour. By clustering the elements around the points of maximum gradient magnitude, the contour’s external energy is lower than it would be if the identical contour is re-sampled to the same number of nodes but with equal node spacing. One possible solution to this problem is to integrate the external energy with respect to the normalised path length l(s) given by:

(2.14)

where L is the total contour length

Note: This sets the incremental path length dl/ds to be scale-independent. If the division by L in (2.14) is omitted then it has been found that the contour becomes unstable and may generate ‘pseudopodia’ especially where there are gently sloping regions of the external energy function.

The normalisation suggested by Cohen (1991) and stated previously in §2.2.8 should not suffer from the problem of collapsing elements. Karaolani et al. (1990) propose a penalty term to attempt to balance the forces causing this effect. Theoretically, the elastic energy term prefers elements to be at least equally spaced (the energy is a sum of squares of node-lengths) and the bending energy term promotes straight and equal-length elements. However, increasing the internal parameters also has the effect of causing the contour to shrink and even collapse totally.

 

2.2.11 Modified Internal Energy and Forces

Another disadvantage of the original active contour model given by (2.1) is that both of the internal energy terms, i.e. the tension and stiffness terms, are scale-variant so that the internal energy can always be reduced by a contraction of the contour. With reference to closed contours, Xu et al. (1994) propose that the internal force components normal to the contour that cause shrinkage be counteracted by adding a weighted internal pressure term similar to that previously used by Sullivan et al. (1990) and Cohen (1991). Sullivan reported the use of variable tension coefficients so that under the action of internal pressure the default shape of the active contour was elliptical. Cohen proposed the pressure term as an active component that was used to inflate the contour from a small trivial initialisation. A problem with the proposal of Xu is that if one requires control of stiffness in a contour, then it is necessary to allow production of bending forces normal to the boundary, which straighten or at least smooth parts of the contour. At certain points along the contour these forces will be generally oriented into the object and at other points the forces will be directed out from the object. However, the notion of using an internal varying pressure-like force to prevent collapse of a closed contour is a useful addition to the model. One may use the internal pressure term to balance only the contraction effect of the tension forces, ignoring any contraction effect of the stiffness forces. The default shape defined by a pressurised contour with non-zero internal energy is a circle. Sullivan et al. (1990) state that this makes the contour behave like a bubble where the pressure to tension ratio is proportional to the area. It is thus possible to calculate the internal pressure required that would maintain contour equilibrium at a pre-defined area. A third method proposed by Radeva and Marti (1995) is to base the internal energy terms on the difference between a static reference contour and the active contour. The following internal energy terms replace the ones given in (2.1):

(2.15)

where v0 is the static reference contour.

The advantage of this method over the pressurised ‘balloon’ or ‘bubble’ model is that the contour defaults to the shape given by the reference contour and not a circle. Keeping the reference contour static means that the model will be scale-variant.

 

2.2.12 Region Pressure Energy

Ivins and Porril (1994) reported a region growing version of the active contour model (‘statistical snake’) that replaced the external gradient energy with a pressure energy term weighted locally along the contour by some distance metric of the external potential field at each point. Using a gray level metric the contour model could be converted into a centroid-linkage region growing algorithm with boundary smoothness constraints and this could be further generalised to respond to texture energy potentials (Laws, 1979) by using the Mahalanobis distance as a metric. If the external energy field was defined to be constant everywhere, the region energy pressure simply equalled the classical pressure energy. Thus the ‘statistical snake’ can be regarded as a generalisation of the pressurised snake. Whilst this is a powerful development of the capability of active contour models it clearly suffers from the same problems that beset standard region growing algorithms when applied to wounds: the wound region is required to have uniform gray-level, colour or textural properties and this is not generally the case.

 

2.2.13 Directional External Energy and Forces

The external energy term defined in (2.1) treats the potential function in an isotropic manner. It is sometimes undesirable to allow all edges to be able to attract a neighbouring part of the contour since they may represent the boundary of an object separate to that being delineated. This is the case in many wound images where strong edges are created by the contrast between the edges of a limb and the background of the image. If part of the contour is placed too close to such an edge it will be deformed to fit the edge of the limb rather than that of the wound. Several authors have formulated directional energy terms to enable the contour model to discriminate between edges based upon their orientation relative to the contour. In the case of wound images it is often found that the wound is darker than the surrounding skin and the surrounding skin lighter than the background, so that gradients of the potential function directed normally outwards from the contour are considered valid wound edges and gradients oriented contrariwise are considered spurious. Radeva and Marti (1995) use distance propagation to create the potential function, so that the potential at any point in the image is a measure of the Euclidean distance to the nearest detected edge pixel. This distance measure is given a polarity depending upon the orientation of the nearest detected edge point. They derive forces from this potential function, given by:

(2.16)

where f is the angle of the vector normal to contour, directed into the (closed) contour and Px,Py are the spatial derivatives of P.

This formulation provides the contour model with the ability to be repelled away from gradients with incorrect orientation, although it does not appear to be rotationally symmetric. The function may be modified to allow a degree of intransigence towards such gradients instead of the repulsion. A possible disadvantage of (2.16) is that it is not expressed in terms of energies but in terms of modified force vectors. There does not seem to be a tangible potential function whose partial derivatives yield the force vector function given, although this should not be taken to imply that such a function does not exist. The lack of an energy formulation reduces the number of options available for conducting the energy minimisation process; i.e. one cannot use a greedy algorithm (§2.2.9) or the Rayleigh-Ritz finite element method.

A directionally weighted potential energy function is reported by Lai and Chin (1993). Since they use a greedy algorithm, which directly minimises energy by local searches, only the energy function needs to be specified. Their potential energy function is given by the following:

(2.17)

where ji is the angle between the image gradient vector at snaxel vi and the corresponding tangent vector to contour, and ½ÑI(vi)½ is the magnitude of the image gradient normalised to [0,1]. The authors do not state whether this normalisation relates to the whole image or just to the points with which the contour is in contact at each iteration step.

 

2.2.14 Setting the Contour Regularisation Parameters

It was stated in §2.2.4 that the setting of the contour regularisation parameters, a and b, which control the elastic and bending energy terms of the contour respectively is problematic. Terzopoulos (1987) has suggested that the parameters be dynamically adjusted along the contour at each iteration in order to conform the contour towards a prescribed ‘natural’ state having a ‘natural arc length’, L(s) and a ‘natural curvature’, C(s). The formulation is given in the continuous case and the parameter (functions) are set by the following:

(2.18)

(2.19)

where |v(s)| is the incremental arc length of the contour and k(s) is the curvature of the contour.

Leymarie and Levine (1993) state that altering the parameters at each iteration using (2.18) and (2.19) has the tendency to make the snaxels oscillate around the prescribed values of natural arc length and natural curvature. This has clear implications for the stability and convergence of the model. They propose a clipping function to constrain the values of a to remain within a prescribed band of acceptable values. The curvature term poses a much more difficult problem (Leymarie and Levine, 1989). To avoid this problem they suggest setting b to a small positive constant, giving the contour a weak tendency to straighten.

Williams and Shah (1992) compare several discrete measures of curvature and propose the following curvature estimate as the most representative value of discrete curvature:

(2.20)

where ui=vi-vi-1 is a vector joining two neighbouring snaxels and q is the angle between two such vectors sharing a common snaxel.

They propose setting a=1 and b=1, with the gradient term weighted by g=1.2 (the gradient energy is normalised locally so that it cannot exceed 1). The stiffness parameter is set to zero at points where the discrete curvature is a local maximum along the contour, provided that (a) its value is greater than a minimum curvature threshold, and (b) the edge strength at that point is above a gradient threshold. This allows the contour to develop a corner at such a point.

Lai and Chin (1993) make use of the minimax principle (Gennert and Yuille, 1988), which they show obviates the need to explicitly set the parameters. The application of the minimax principle assumes that the energy of an active contour model in equilibrium is a continuous convex function of the regularisation parameters. It proposes that the maximum point on the surface of this function corresponds to the parameter levels that provide the best trade-off between noise-insensitivity (high regularisation) and the need to allow enough flexibility of the contour to adequately follow the desired contour (low regularisation). Applying the minimax principle, the energy term to be minimised at each iteration is the maximum of the energy terms that constitute the model.

 

2.2.15 Initialisation of Active Contour Models

The iterative equations for solving an active contour problem require an initial approximation to the desired solution. In this work, the initialisation relies upon a manual delineation of the boundary so that the active contour model is used to refine the solution proposed by a trained delineator.

Cohen (1991) describes a ‘balloon’ active contour model that uses an internal pressure force – which by definition acts with constant pressure normal to the boundary everywhere – to inflate a contour situated inside the boundary of the desired object. The pressure and image potential forces were controlled so that noise pixels did not trap the contour before the desired boundary was located. The need to simply initialise the contour to be a small region within the region of interest is an attractive property of this model. Unfortunately, in the context of wound measurement, the model possesses two overriding disadvantages:

Lai and Chin (1995) show that the Generalised Hough Transform (GHT) can be used to model a contour. They use the GHT to obtain an approximate solution and state that the GHT essentially applies a rigid template model to the problem. The active contour model (deformable template) then refines the result of the GHT. This approach is most useful when the object(s) of interest are of similar shape such that they can be adequately modelled – if only approximately – by the GHT. However, in the case of wounds, the shapes are too diverse for a GHT to be applicable. However, it may be that once a wound has been delineated and this delineation refined by application of an active contour model, then the shape information gained by so doing could be employed to refine subsequent wound measurements and even to subsequently dispense with the need for delineation. If an initial wound template could be used as a reliable deformable model for subsequent measurements then the g-snake technique (Lai and Chin, 1995) could be used to provide more robust and automated wound measurement. To implement this requires a controlled study of a large enough set of wounds to ensure the broadness of its applicability.

The goal of automated initialisation could be pursued by using region-based methods such as pixel classification schemes based on colour properties of the wound components and surrounding skin. As previously stated, these properties vary from one patient to another so that supervised classification is not flexible enough to allow discrimination in many cases. Clustering techniques are better suited to the problem, although it is foreseeable that there will be very many instances where such initialisation will not place the entire contour in the vicinity of the wound’s boundary, since there are ‘false’ edges outside and more commonly inside the wound.

 

2.2.16 Global Energy Minimum Seeking

An active contour model algorithm attempts to iterate to convergence at a local energy minimum. Regularisation of the contour itself and also of the potential image function suppresses the many noisy energy minima that exist in the vicinity of each true solution. In an image with other edges beside those belonging to the object of interest, other local energy minima will exist and if the contour model is not correctly seeded, will converge to the wrong solution and commonly, to a hybrid solution with parts of various different boundaries forming the whole.

If enough information can be provided to make the energy of the true wound contour the global minimum out of all the energy minima existing in the domain of the problem, one may attempt to seek such a minimum using a trivial initialisation, e.g. a circle centred at the centre of the image of some default size, perhaps relative to the image itself. This implies that not only is the energy model both specific to wounds and general enough to encompass all wounds but also that the algorithm used to fit the model can find the global minimum. One method of finding a global minimum is that of simulated annealing (Aarts and Laarhoven, 1987). Some experimentation with this method has shown that a major factor in the likely success of this approach is an appropriate parameterisation of the contour. Storvik (1992) used simulated annealing to solve an active contour optimisation problem, but allowed only localised deformations to be made to the contour.

Staib and Duncan (1992) parameterised the boundary of an object using a Fourier decomposition which generalised to be scale and rotation invariant. They used an independent multi-dimensional Gaussian probability distribution to model the expected variation in the Fourier-decomposed parameter vector elements for their search shape. Using Bayes’ theorem the boundary finding model is cast as a maximum a-posteriori problem. Taking the maximum a-priori boundary vector as an initial solution they used a gradient-descent algorithm to converge to, firstly, a coarse-scale solution achieved by blurring the image with a Gaussian filter, then using the converged result as a re-initialisation, a second finer-scale solution was sought. Use of the Gaussian probability function implies that the mean and variance of the parameters for a particular domain of possible templates or corresponding object boundaries are known.

In the case of finding a wound’s boundary, the parameter distributions could be estimated from the manually delineated boundaries of a set of wounds. This process is likely to require many training samples of wound boundaries, especially if the number of parameters required to model the distributions is large. Specifying a distribution of the parameters covering all encounterable wounds is likely to yield a vague description of any one wound case.

 

2.3 Summary

This chapter has appraised the various segmentation algorithms that are commonly reported in the available literature. This work maintains that the varied visual properties of a substantial proportion of leg ulcer images limits the applicability of many algorithms which require the data representing each image region to be homogeneous with respect to some criterion. Indeed the principle of dividing an image into homogeneous regions is commonly cited by standard texts on image processing. Leg ulcers tend to be somewhat varied and slough formed from dead tissue complicates their appearance. In addition to this, epithelialisation tissue and the quality of the surrounding ‘undamaged’ skin is varied among different wounds. However, the leg ulcer images have been gathered from several different sources and under varied illumination conditions. Improving the stability and homogeneity of illumination is certainly a worthwhile avenue of further work. The active contour models appear to be the most widely applicable technique for measuring the size of a leg ulcer. There is no assumption that the wound is homogeneous; instead the requirements for successful measurement are fairly good contrast around the wound (especially at points of high boundary curvature) and good initialisation. The following chapter formulates four active contour model algorithms that incorporate many of the proposals made by other researchers to improve, adapt and extend the capabilities of active contour models.


Title Page

4 Manual Delineation

Appendix A

Contents

5 Parameter Setting

Appendix B

1 Introduction

6 Results

Appendix C

2 Literature Survey

7 Discussion

Colour Plates

3 Algorithm Development

8 Conclusion

References